##################################################################
#
# REJECTION METHOD = CLASSROOM
#
# Target   : pi(theta) = N(theta,theta0,C0)*N(x,theta,1)
# Proposal : q(theta) = t(x,5,10)
#
##################################################################
#
# 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/
#
#################################################################
pi = function(theta){exp(-0.5*(x-theta)^2)*exp(-0.5*(theta-theta0)^2/C0)}
q  = function(theta){dt((theta-x)/5,10)}

# Prior hyperparameters
theta0 = 0
C0     = 1

# Observed data
x      = 1

# What is a possible value of A, the envelop constant?
As = c(1,2,3)
thetas = seq(-30,30,length=1000)
pis = pi(thetas)
qs  = q(thetas)
par(mfrow=c(2,3))
for (A in As){
  L = min(pis,A*qs)
  U = max(pis,A*qs)
  plot(thetas,pis,type="l",xlab="",ylab="Density",main="",ylim=c(L,U))
  lines(thetas,A*qs,col=2,lwd=3)
  title(paste("A=",A,sep=""))
}
for (A in As){
  plot(thetas,pis/(A*qs),type="l",xlab="",ylab="pi/A*q",main="")
  abline(h=1,col=2,lwd=3)
}

# Rejection method in action
# --------------------------
A       = 2.0                        # Setting envelop constant
M       = 10000                      # Monte Carlo size
thetas  = x+sqrt(5)*rt(M,10)         # Drawing from proposal q(.) 
prob    = pi(thetas)/(A*q(thetas))   # Computing acceptance probabilities
u       = runif(M)                   # Drawing the uniforms
thetas1 = thetas[u<prob]             # Draws from pi(.) (the target distribution)

par(mfrow=c(1,1))
breaks = seq(min(thetas1),max(thetas1),length=30)
hist(thetas1,prob=TRUE,main="",breaks=breaks)
ths = seq(x/2-3*sqrt(1/2),x/2+3*sqrt(1/2),length=1000)
lines(ths,dnorm(ths,x/2,sqrt(1/2)),col=2,lwd=2)
title(paste("Rejection = ",round(100-100*length(thetas1)/M,1),"%",sep=""))