#############################################################################################################
#
# Example 7.10 (pages 274-5) – Figure 7.5 (page 275)  
#
# Random walk Metropolis and delayed rejection Metropolis algorithms to sample from 
#
#                  pi(theta)  =  0.9f_N(theta;0,1)  +  0.1f_N(theta;10,1)
#
##############################################################################################################
#
# Author : Hedibert Freitas Lopes
#          Graduate School of Business
#          University of Chicago
#          5807 South Woodlawn Avenue
#          Chicago, Illinois, 60637
#          Email : hlopes@ChicagoGSB.edu
#
###############################################################################################################
logpost = function(theta){log(0.9*dnorm(theta)+0.1*dnorm(theta,10,1))}

# Common settings
# ---------------
burnin  = 0
M       = 15000
Vtheta  = 2
theta0  = 10.5

# Random Walk Metropolis-Hastings
# -------------------------------
set.seed(524157)
samples = NULL
theta   = theta0
sample  = NULL
for (iter in 1:(burnin+M)){
  th1  = rnorm(1,theta,Vtheta)
  prob = exp(logpost(th1)-logpost(theta))
  if (runif(1) < prob){
    theta = th1
  }
  sample = c(sample,theta)
}
sample = sample[(burnin+1):(burnin+M)]

# Delayed Rejection Metropolis Algorithm 
# --------------------------------------
Vtheta1 = 2*Vtheta
theta   = theta0
sample1 = NULL
for (iter in 1:(burnin+M)){
  th1   = rnorm(1,theta,Vtheta)
  prob  = min(1,exp(logpost(th1)-logpost(theta)))
  u     = runif(1)
  ind   = u < prob
  if (ind){
    theta = th1
  }else{
    th2   = rnorm(1,th1,Vtheta1)
    prob1 = min(1,exp(logpost(th1)-logpost(th2)))
    num   = logpost(th2)+dnorm(th1,th2,Vtheta,log=TRUE)+log(1-prob1)+dnorm(theta,th1,Vtheta1)
    den   = logpost(theta)+dnorm(th1,theta,Vtheta,log=TRUE)+log(1-prob)+dnorm(th2,th1,Vtheta1)
    u     = runif(1)
    ind   = u < exp(num-den)
    if (ind){
      theta = th2
    }  
  }
  sample1 = c(sample1,theta)
}
sample1 = sample1[(burnin+1):(burnin+M)]

# Figure  7.5
# -----------
theta = seq(-4.1,15,length=1000)
ftheta = exp(logpost(theta))
breaks = seq(-4.1,15,length=50)
par(mfrow=c(3,2))
plot(seq(1,burnin+M,by=20),sample[seq(1,burnin+M,by=20)],xlab="iterations",ylab="",main="Random walk Metropolis",ylim=range(c(sample,sample1)),axes=F,type="l")
axis(1,at=seq(1000,burnin+M,by=2000),label=seq(1000,burnin+M,by=2000))
axis(2,at=seq(-3.5,10.5,length=5),label=seq(-3.5,10.5,length=5))
plot(seq(1,burnin+M,by=20),sample1[seq(1,burnin+M,by=20)],xlab="iterations",ylab="",main="Delayed rejection Metropolis",ylim=range(c(sample,sample1)),axes=F,type="l")
axis(1,at=seq(1000,burnin+M,by=2000),label=seq(1000,burnin+M,by=2000))
axis(2,at=seq(-3.5,10.5,length=5),label=seq(-3.5,10.5,length=5))
acf(sample,xlab="lag",ylab="",main="",axes=F)
axis(1)
axis(2)
acf(sample1,xlab="lag",ylab="",main="",axes=F)
axis(1)
axis(2)
hist(sample,breaks=breaks,xlab="",ylab="",main="",prob=T,axes=F,ylim=c(0,0.4))
axis(2)
axis(1,at=seq(-3.5,10.5,length=5),label=seq(-3.5,10.5,length=5))
lines(theta,exp(logpost(theta)))
hist(sample1,breaks=breaks,xlab="",ylab="",main="",prob=T,axes=F,ylim=c(0,0.4))
axis(2)
axis(1,at=seq(-3.5,10.5,length=5),label=seq(-3.5,10.5,length=5))
lines(theta,exp(logpost(theta)))