# Algo: epage 70
# Problem: epage 75
#
# 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):
#
M = sqrt( 2 * pi ) / sqrt(exp(1))

Nsim = 2500

fFunc = function(x){
  dnorm(x)
}

x=NULL
while( length(x)<Nsim ){
  # generate samples from the candidate density
  y = rcauchy( Nsim*M, scale=1. )
  x = c(x, y[ runif(Nsim*M)*M < fFunc(y)/dcauchy(y) ])
}
x = x[1:Nsim]

hist(x,freq=F)
xRng = seq( min(X), max(X), length.out=100 )
lines(xRng,fFunc(xRng),lwd=2,col="sienna")



# Part (b):
#

# plot the two functions before any scaling:
# 
x = seq(0,100,length.out=100)
plot( x, dgamma(x,4.3,1/6.2) ) # this if f(x) the density we want 
lines( x, dgamma(x,4,1/7) ) # this is g(x) the density we have 

M = optimize(f=function(x){
  num = dgamma(x,4.3,1/6.2)
  den = dgamma(x,4,1/7)
  num/den 
}, interval=c(0,+20.),maximum=T)$objective

# plot the two functions after scaling:
# 
postscript("../../WriteUp/Graphics/Chapter2/ex_2_20_f_N_Mg.eps", onefile=FALSE, horizontal=FALSE)

plot( x, dgamma(x,4.3,1/6.2) )
lines( x, M*dgamma(x,4,1/7) )

dev.off()

X=NULL
while( length(X)<Nsim ){
  Y = rgamma( M*Nsim, 4,1/7 )
  X = c(X, Y[ M*runif(M*Nsim) < dgamma( Y,4.3,1/6.2 )/dgamma( Y, 4,1/7 ) ] )
}

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

hist(X,freq=F)
xRng = seq( min(X), max(X), length.out=100 )
lines(xRng,dgamma(xRng,4.3,1/6.2),lwd=2,col="sienna")

dev.off()