#
# 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.
# 
#-----

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

Nsims = 10^4
samps = runif(Nsims,min=0,max=1./20) # random draws

integrand = (1/sqrt(2*pi)) * exp( -1/(2*samps^2) )/(samps^2)
I_est = mean( integrand )
print( sprintf("Problem 3.3: N= %8d; MC estimate= %10.6g", Nsims, I_est) )

# To study convergence:
#
integrand_fn = function(x) { (1/sqrt(2*pi)) * exp( -1/(2*x^2) )/(x^2) }

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

par( mfrow=c(2,1) )
curve( integrand_fn, from=0, to=1./20, xlab="independent variable (v)", ylab="integrand", lwd=2 )
print( integrate( integrand_fn, 0, 1./20 ) )

x = integrand_fn( runif(Nsims)/20 )
estint = cumsum(x)/1:Nsims
esterr = sqrt( cumsum((x - estint)^2) )/(1:Nsims)
plot( estint, type="l", xlab="number of samples", ylab="mean and error range", lwd=2, ylim=mean(x)+20*c( -esterr[Nsims], +esterr[Nsims] ) )
lines( estint - 2 * esterr, col="gold", lwd=2 )
lines( estint + 2 * esterr, col="gold", lwd=2 )

par( mfrow=c(1,1) )

#dev.off()