#Target distribution 
target = function(theta){
  0.2*dnorm(theta,2,1)+0.8*dnorm(theta,10,1)
}
# 1st Proposal distributon
proposal1 = function(theta){
  dnorm(theta,6,5)
}

# 2nd Proposal distributon
proposal2 = function(theta,thetaold){
  dnorm(theta,thetaold,2)
}

theta = seq(-9,21,length=1000)
par(mfrow=c(1,1))
plot(theta,target(theta),xlab=expression(theta),ylab="Density",type="l",lwd=2)
lines(theta,proposal1(theta),lwd=2,col=2)
legend("topleft",legend=c("Target","Proposal"),col=1:2,lwd=2,lty=1)

# Metropolis-Hastings algorithm (based on 1st proposal)
th     = -9
set.seed(233243)
M0     = 5000
M      = 5000
niter  = M0+M
thetas = rep(0,niter)
for (i in 1:niter){
  thstar = rnorm(1,6,5)
  num = target(thstar)/proposal1(thstar)
  den = target(th)/proposal1(th)
  accept = min(1,num/den)
  if (runif(1)<accept)
    th = thstar
  thetas[i] = th
}

# MCMC output and approximation
par(mfrow=c(2,3))
plot(thetas,xlab="iterations",ylab="draws",main="Metropolis-Hastings draws\n Independent proposal")
abline(v=M0,col=2)
acf(thetas[(M0+1):niter],main="")
breaks=seq(min(thetas),max(thetas),length=30)
hist(thetas[(M0+1):niter],prob=TRUE,xlab=expression(theta),main="",xlim=range(theta),breaks=breaks)
lines(theta,target(theta),lwd=2)
box()
title("Target approximation")



# Metropolis-Hastings algorithm (based on 1st proposal)
th     = -9
set.seed(233243)
M0     = 5000
M      = 5000
niter  = M0+M
thetas = rep(0,niter)
for (i in 1:niter){
  thstar = rnorm(1,th,2)
  num = target(thstar)/proposal2(thstar,th)
  den = target(th)/proposal2(th,thstar)
  accept = min(1,num/den)
  if (runif(1)<accept)
    th = thstar
  thetas[i] = th
}

# MCMC output and approximation
plot(thetas,xlab="iterations",ylab="draws",main="Metropolis-Hastings draws\n Random walk proposal")
abline(v=M0,col=2)
acf(thetas[(M0+1):niter],main="")
breaks=seq(min(thetas),max(thetas),length=30)
hist(thetas[(M0+1):niter],prob=TRUE,xlab=expression(theta),main="",xlim=range(theta),breaks=breaks)
lines(theta,target(theta),lwd=2)
box()
title("Target approximation")