#
#
# Written by:
# -- 
# John L. Weatherwax                2009-04-21
# 
# email: wax@alum.mit.edu
# 
# Please send comments and especially bug reports to the
# above email address.
# 
#-----

# Part (a):
# 

# Find the largest value of f(y)/g(y) (note the Gaussian exponential factor cancels):
# 
M = optimize(f=function(x){
  num = ( sin(6*x)^2 + 3*( cos(x)^2 ) * (sin(4*x))^2 + 1 )
  den = ( 1/sqrt(2*pi) )
  num/den 
}, interval=c(0,+1.),maximum=T)$objective


# Here we will plot the density f and the upper bounding density $M g(x)$.
# 
###postscript("../../WriteUp/Graphics/Chapter2/ex_2_18_f_function.eps", onefile=FALSE, horizontal=FALSE)
x = seq(-5,+5,length.out=1000)
f = exp(-0.5*x^2)*( sin(6*x)^2 + 3*( cos(x)^2 ) * (sin(4*x))^2 + 1 )

plot( x, f )
lines( x, f, lwd=2, col="black")

g = M * exp(-0.5*x^2)/sqrt(2*pi)
lines( x, g )

###dev.off() 

# Generate Nsim random variables from f and estimate the acceptance rate of this proceedure 
#
Nsim = 2500

fFunc = function(x){
  exp(-.5*x^2)*( sin(6*x)^2 + 3*( cos(x)^2 ) * (sin(4*x))^2 + 1 )
}
gFunc = function(x){
  exp(-.5*x^2)/sqrt( 2*pi ) 
}

X = NULL
U = NULL
numDrawsFound = 0
while( numDrawsFound<Nsim ){
  X = c(X,rnorm(Nsim)) # extract some draws from our proposal distribution 
  U = c(U,runif(Nsim)) # extract the corresponding uniform random variable 
  numDrawsFound = sum( fFunc(X)/(M*gFunc(X)) >= U ) # find out which draws passed the accept-reject threshold
}

prob_acceptance = numDrawsFound/length(X)
hat_f = 1/(prob_acceptance*M) # the nomalizing factor for the f(x) density 

# extract out the true samples from f(x) and estimate the probability of acceptance (1/M):
#
inds = fFunc(X)/(M*gFunc(X)) >= U 
X = X[inds]

###postscript("../../WriteUp/Graphics/Chapter2/ex_2_18_normalized_f_density.eps", onefile=FALSE, horizontal=FALSE)

hist(X,freq=F,breaks=100)
xRng = seq( min(X), max(X), length.out=100 )
#lines(xRng,fFunc(xRng)/(prob_acceptance*M),lwd=2,col="sienna")
lines(xRng,hat_f * fFunc(xRng),lwd=2,col="sienna")

###dev.off()