#############################################################################################################
#
#  Example 1.3 (page 27) 
#  Figure 1.3 (page 27) 
#
#  Rejection method: sampling N(0,1) draws from Cauchy(0,2.5) draws. 
#
##############################################################################################################
#
# Author : Hedibert Freitas Lopes
#          Graduate School of Business
#          University of Chicago
#          5807 South Woodlawn Avenue
#          Chicago, Illinois, 60637
#          Email : hlopes@ChicagoGSB.edu
#
###############################################################################################################
dtv = function(x,df=1,mu=0,sig2=1){
  dt((x-mu)/sqrt(sig2),df)/sqrt(sig2)
}
# Figure 1.3
# ----------
set.seed(28386)
x0     = seq(-6,6,length=1000)
f1     = dnorm(x0)
f2     = dtv(x0,1,0,2.5)
A      = sqrt(pi*1.25)
M      = 5000
x      = sqrt(2.5)*rt(M,1)
weight = dnorm(x)/(A*dtv(x,1,0,2.5))
ind    = runif(M)<weight
x1     = x[ind]
x1     = sort(x1)

par(mfrow=c(1,2))
plot(x0,f1/(A*f2),xlab="(a)",ylab="",main="",type="l",axes=F);axis(1);axis(2)
hist(x1,xlab="(b)",ylab="",prob=T,breaks=seq(min(x1),max(x1),length=25),main="")
lines(x0,f1)