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

Nsim = 10^6;

alpha = 1

M = sqrt( 2/pi ) * exp(1/2)

randg = function(Nsim,alpha){
  # a function to generate Nsim random variable from a double exponential using the inverse transform 
  # 
  X = double(Nsim) 
  U = runif(Nsim)
  X[ U <= 1/2 ] =  log(2*   U[ U <= 1/2 ])/alpha 
  X[ U  > 1/2 ] = -log(2*(1-U[ U >  1/2 ]))/alpha
  X
}


# draw samples from g(y) using the known probability of acceptance
#
x = NULL
while( length(x)<Nsim ){
  Y = randg(Nsim*M,alpha)
  
  fY = dnorm( Y, mean=0, sd=1 )
  gY = (alpha/2) * exp(-alpha*abs(Y))
  r  = fY/(M*gY) 

  x = c(x,Y[runif(Nsim*M)<=r])
}
x = x[1:Nsim] # these should be Gaussian N(1,1)

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

par(mfrow=c(1,2))
Z = rnorm(Nsim)
zmin = min(min(Z),min(x))
zmax = max(max(Z),max(x))
b = seq(zmin,zmax,length.out=100)
hist(Z,freq=F,main="Using rnorm",breaks=b)
hist(x,freq=F,main="Using Double Exponential",breaks=b)

##dev.off()