#
#
# 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 here
  (1/sqrt(2*pi)) * exp( -0.5*(x - mu_bar)^2 )
}
fFunc = function(x){ # don't do the truncation when x<a here
  exp( -0.5*x^2 )
}

Nsim = 2500

x=NULL
while( length(x)<Nsim ){
  
  y = rnorm( Nsim, mean=mu_bar, sd=1. ) # generate samples from the candidate density

  yKeep = y[ y > 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<a here
  alpha * exp( -alpha*(x - a) )
}
fFunc = function(x){ # don't do the truncation when x<a here
  exp( -0.5*x^2 )
}

Nsim = 2500

x=NULL
while( length(x)<Nsim ){
  
  y = rexp( Nsim, rate=alpha ) # generate samples from the candidate density
  y = y + a # now tranlate 

  x = c(x, y[ runif(Nsim)*MTilde < fFunc(y)/gFunc(y,alpha,a) ])
}
x = x[1:Nsim]


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

# Plot this as a normalized density:
# 
hist(x,freq=F,main="truncated normal dist. via translated exponential")
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()