#
# EPage 82
# 
# 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) plot the integrands (numerator and denominator):
#
thetas = seq( -10, +10, length.out=100 )

#postscript("../../WriteUp/Graphics/Chapter3/ex_3_1_numerator_integrand_plots.eps", onefile=FALSE, horizontal=FALSE)
#postscript("../../WriteUp/Graphics/Chapter3/ex_3_1_denominator_integrand_plots.eps", onefile=FALSE, horizontal=FALSE)

par(mfrow=c(1,3))
for( x in c(0,2,4) ){
  integrand = ( thetas/( 1 + thetas^2 ) ) * exp( - ( x - thetas )^2 / 2 )
  plot( thetas, integrand, type="l", main=sprintf("numerator (x=%10.f)",x) )

  #integrand = ( 1/( 1 + thetas^2 ) ) * exp( - ( x - thetas )^2 / 2 )
  #plot( thetas, integrand, type="l", main=sprintf("denominator (x=%10.f)",x) )
}
par(mfrow=c(1,1))

#dev.off()


# Estimate the integrals using MC by drawing samples from a Cauchy distribution: 
# 
set.seed(0)

print( "MC delta(x) estimation using cauchy draws ..." )
for( x in c(0,2,4) ){
  Nsims = 10^4
  r_thetas = rcauchy(Nsims) # random draws
  num = mean( r_thetas * dnorm(r_thetas,mean=x) ) # numerator
  den = mean( dnorm(r_thetas,mean=x) ) # denominator
  print( sprintf("x= %3.f; N= %8d; MC estimate= %10.6f", x, Nsims, num/den) )
}


# Part (b): Estimate the convergence of our ratio by considering the convergence of each factor in the ratio:
#

cauchy_est = function(x,Nsims){
  # Estimates the integral in the numerator of delta(x) using draws from a Cauchy distribution
  #
  r_thetas = rcauchy(Nsims) # random draws
  integrand = r_thetas * dnorm(r_thetas,mean=x) # numerator 
  int_est = cumsum( integrand )/(1:Nsims)
  int_err = sqrt( cumsum( (integrand - int_est)^2 ) )/(1:Nsims)
  
  list( int_est, int_err ) 
}

normal_est = function(x,Nsims){
  # Estimates the integral in the numerator of delta(x) using draws from a normal distribution
  #
  r_thetas = rnorm(Nsims,mean=x) # random draws
  integrand = r_thetas * dcauchy(r_thetas) # numerator 
  int_est = cumsum( integrand )/(1:Nsims)
  int_err = sqrt( cumsum( (integrand - int_est)^2 ) )/(1:Nsims)
  
  list( int_est, int_err ) 
}


#postscript("../../WriteUp/Graphics/Chapter3/ex_3_1_numerator_bayes_convergence_plots.eps", onefile=FALSE, horizontal=FALSE)

par(mfrow=c(1,3))
for( x in c(0,2,4) ){
  Nsims = 10^3
  
  c_expts = cauchy_est(x,Nsims)
  v  = c_expts[[1]]
  se = c_expts[[2]]
  plot( 1:Nsims, v, type="l", col="black", main=sprintf("estimates of numerator (x=%3.f)",x) )
  lines( v + 2*se, col="red" )
  lines( v - 2*se, col="red" )

  n_expts = normal_est(x,Nsims)
  v  = n_expts[[1]]
  se = n_expts[[2]]
  lines( v, type="l", col="gray" )
  lines( v + 2*se, col="pink" )
  lines( v - 2*se, col="pink" )
}
par(mfrow=c(1,1))

#dev.off()


# Part (c): The same experiments but with draws from a random normal rather than the cauchy 
#
set.seed(0)

print( "MC delta(x) estimation random normal draws ..." )
for( x in c(0,2,4) ){
  Nsims = 10^4
  r_thetas = rnorm(Nsims,mean=x) # random draws
  num = mean( r_thetas * dcauchy(r_thetas) )
  den = mean( dcauchy(r_thetas) )
  print( sprintf("x= %3.f; N= %8d; MC estimate= %10.6f", x, Nsims, num/den) )
}