#
#
# 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.
#
#-----
a = 3 # a large value that makes naive simulation from the truncated normal distribution hard
# Part (c):
#
M = sqrt( 2 * pi ) * exp( -0.5*a^2 ) # this is *not* the reciprical of the acceptance rate since the densities are not correctly normalized
mu_bar = a
# Plot the two densities f(x) and M*g(x) to make sure that f(x) < M*g(x):
#
xRange = seq( a, 4*a, length.out=100 )
f = exp( -0.5*xRange^2 )
plot( xRange, f )
g = (1/sqrt(2*pi)) * exp( -0.5*(xRange - mu_bar)^2 )
lines( xRange, M*g )
gFunc = function(x,mu_bar){ # don't do the truncation when x a ] # these are possible candidates (other points have f(y)=0 will not pass the acceptance phase)
nKeep = length(yKeep)
y = yKeep
x = c(x, y[ runif(nKeep)*M < fFunc(y)/gFunc(y,mu_bar) ])
}
x = x[1:Nsim]
##postscript("../../WriteUp/Graphics/Chapter2/ex_2_21_pt_trunc_norm_dist.eps", onefile=FALSE, horizontal=FALSE)
# Plot this as a normalized density:
#
hist(x,freq=F,main="truncated normal dist. via normal")
xRng = seq( min(x), max(x), length.out=100 )
normConstant = sqrt(2*pi) * ( 1 - pnorm(a) )
lines(xRng,fFunc(xRng)/normConstant,lwd=2,col="sienna")
##dev.off()
# Part (d):
#
alpha = ( a + sqrt( a^2 + 4 ) )/2 # the optimal value for the translated exponential density
MTilde = (1/alpha)*exp( 0.5*alpha^2 - alpha*a ) # the optimal value for MTilde
# Plot the two densities f(x) and M*g(x) to make sure that f(x) < M*g(x):
#
xRange = seq( a, 4*a, length.out=100 )
f = exp( -0.5*xRange^2 )
plot( xRange, f )
g = alpha * exp( -alpha*(xRange - a) )
lines( xRange, M*g )
gFunc = function(x,alpha,a){ # don't do the truncation when x