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

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,proposal(theta),lwd=2,col=2)
legend("topleft",legend=c("Target","Proposal"),col=1:2,lwd=2,lty=1)

# Initial value
th     = -9

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

# MCMC output and approximation
par(mfrow=c(2,2))
plot(theta,target(theta),xlab=expression(theta),ylab="Density",type="l",lwd=2)
lines(theta,proposal(theta),lwd=2,col=2)
legend("topleft",legend=c("Target","Proposal"),col=1:2,lwd=2,lty=1,bty="n")
plot(thetas,xlab="iterations",ylab="draws",main="Metropolis-Hastings draws")
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")