######################################################
#
#  Unit 2 - Example iii.
#
#  Rejection method
#
#      Sampling N(0,1) draws from Cauchy(0,2.5) draws. 
#      Sampling N(0,1) draws from U(-6,6) 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)
}

M      = 2000
set.seed(28386)
par(mfrow=c(2,2))

# Cauchy proposal
# ---------------
x0     = seq(-6,6,length=1000)
f1     = dnorm(x0)
f2     = dtv(x0,1,0,2.5)
A      = sqrt(pi*1.25)
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)
plot(x0,f1/(A*f2),xlab="PROPOSAL: CAUCHY(0,2.5)",ylab="",main=expression(pi(theta)/q(theta)),type="l",axes=F);axis(1);axis(2)
hist(x1,xlab=paste("Acceptance rate=",round(100*length(x1)/M,2),"%",sep=""),ylab="",prob=T,breaks=seq(min(x1),max(x1),length=25),main="")
title(paste("A=",round(A,3),sep=""))
lines(x0,f1)

# Uniform proposal
# ----------------
x0     = seq(-6,6,length=1000)
f1     = dnorm(x0)
f2     = 1/12
A      = 12/sqrt(2*pi)
x      = -6+12*runif(M)
weight = dnorm(x)*sqrt(2*pi)
ind    = runif(M)<weight
x1     = x[ind]
x1     = sort(x1)

plot(x0,f1/(A*f2),main=expression(pi(theta)/q(theta)),ylab="",xlab="PROPOSAL: UNIFORM(-6,6)",type="l",axes=F);axis(1);axis(2)
hist(x1,xlab=paste("Acceptance rate=",round(100*length(x1)/M,2),"%",sep=""),ylab="",prob=T,breaks=seq(min(x1),max(x1),length=25),main="")
title(paste("A=",round(A,3),sep=""))
lines(x0,f1)