pmix = function(theta){
  pi*dnorm(theta,mu[1],1)+(1-pi)*dnorm(theta,mu[2],1)
}

pdf(file="comparingMCMCschemes.pdf",width=15,height=10)

for (dist in c(0,5,10,20)){
pi     = 0.8
mu     = c(0,dist)
thetas = seq(-5,dist+5,length=1000)
den = pmix(thetas)

par(mfrow=c(1,1),c(1.02,0.82,0.82,0.42))
plot(thetas,den,type="l",xlab=expression(theta),ylab="Density")

# MCMC chain
M0    = 1000
M     = 2000
niter = M0 + M

# Random-walk Metropolis-Hastings
set.seed(3141593)
sds.rwmh   = c(1,5,10,20)
draws.rwmh = matrix(0,niter,4)
par(mfrow=c(2,6),mai=c(0.4,0.4,0.8,0.2))
for (i in 1:4){
  sd.rwmh   = sds.rwmh[i]
  theta     = 5
  draw.rwmh = rep(0,niter)
  for (iter in 1:niter){
    theta1  = rnorm(1,theta,sd.rwmh)
    accept  = min(1,pmix(theta1)/pmix(theta))
    u       = runif(1)
    if (u<accept){
      theta = theta1
    }
    draw.rwmh[iter] = theta
  }
  draws.rwmh[,i] = draw.rwmh
  ts.plot(draw.rwmh,xlab="Iterations",ylab="Draws",main="")
  abline(v=M0,col=2,lwd=2)
  acf(draw.rwmh[(M0+1):niter],main="")
  legend("topright",legend=paste("RW-MH: stdev=",sd.rwmh,sep=""),bty="n")
  hist(draw.rwmh[(M0+1):niter],prob=TRUE,main="",xlim=range(thetas),
       xlab=expression(theta),ylab="Density",
       ylim=range(den))
  lines(thetas,den,col=2,lwd=2)
}

# Independent Metropolis-Hastings
set.seed(2718282)
mean.imh  = 5
sds.imh   = c(0.1,0.5,1,5)
draws.imh = matrix(0,niter,4)
par(mfrow=c(2,6),mai=c(0.4,0.4,0.8,0.2))
for (i in 1:4){
  sd.imh   = sds.imh[i]
  theta    = 5
  draw.imh = rep(0,niter)
  for (iter in 1:niter){
    theta1  = rnorm(1,mean.imh,sd.rwmh)
    accept  = min(1,pmix(theta1)/pmix(theta)*
                    dnorm(theta,mean.imh,sd.rwmh)/
                    dnorm(theta1,mean.imh,sd.rwmh))
    u       = runif(1)
    if (u<accept){
      theta = theta1
    }
    draw.imh[iter] = theta
  }
  draws.imh[,i] = draw.imh
  ts.plot(draw.imh,xlab="Iterations",ylab="Draws")
  abline(v=M0,col=2,lwd=2)
  acf(draw.imh[(M0+1):niter],main="")
  legend("topright",legend=paste("I-MH: stdev=",sd.imh,sep=""),bty="n")
  hist(draw.imh[(M0+1):niter],prob=TRUE,main="",xlim=range(thetas),
       xlab=expression(theta),ylab="Density",ylim=range(den))
  lines(thetas,den,col=2,lwd=2)
}


}

dev.off()