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

# fix the random.seed, so I'll always get the same answer:
set.seed(10131985)

Nsim = 2*10^3 # must be an even number for the Box-Muller transofrm 

# The CLT random normal generator:
# 
U = runif(12*Nsim,min=-0.5,max=0.5)
U = matrix(U,nrow=Nsim,ncol=12)
Z_CLT = apply(U,1,sum)

# The Box-Muller algorithm:
# 
U1 = runif(Nsim/2)
U2 = runif(Nsim/2)

X1 = sqrt( -2*log(U1) ) * cos( 2 * pi * U2 )
X2 = sqrt( -2*log(U1) ) * sin( 2 * pi * U2 )
Z_BM = c(X1,X2)

# The R generator:
#
Z_R = rnorm(Nsim)

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

par(mfrow=c(1,3))
hist(Z_CLT,freq=F,main="Central Limit Theorem",breaks=100)
hist(Z_BM,freq=F,main="Box Muller",breaks=100)
hist(Z_R,freq=F,main="R(rnorm)",breaks=100)

###dev.off()

#plot( Z_CLT, Z_BM )