################################################################################
#
#  Rejection method: sampling N(0,1) draws from Cauchy(0,2.5) draws. 
#                    sampling N(0,1) draws from U(-6,6) draws.
#
################################################################################
#
# HEDIBERT FREITAS LOPES
# Associate Professor of Econometrics and Statistics
# The University of Chicago Booth School of Business
# 5807 South Woodlawn Avenue
# Chicago, Illinois, 60637
# Email : hlopes@ChicagoGSB.edu
# URL: http://faculty.chicagobooth.edu/hedibert.lopes/research/
#
################################################################################

dtv = function(x,df=1,mu=0,sig2=1){
  dt((x-mu)/sqrt(sig2),df)/sqrt(sig2)
}

pdf(file="rejection.pdf",height=6,width=12)
M      = 2000
set.seed(28386)
par(mfrow=c(1,4))

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

dev.off()