# Target   = t_4(0,1)
# Proposal = N(theta(i-1),sd)
set.seed(12345)
M0     = 10000
M      = 10000
niter  = M0+M
thetas = rep(0,niter)

par(mfrow=c(3,3))
for (sd in c(0.1,1,2)){
  theta = 2
  for (iter in 1:niter){
    theta.prop = rnorm(1,theta,sd)
    logA = dt(theta.prop,df=4,log=TRUE)-dt(theta,df=4,log=TRUE)
    logacceptance = min(0,logA)
    if (log(runif(1))<logacceptance){
      theta = theta.prop
    }
    thetas[iter] = theta
  }
  ts.plot(thetas,xlab="Iterations",ylab=expression(theta))
  title(paste("SD=",sd,sep=""))
  abline(v=M0,col=2)
  acf(thetas[(M0+1):niter],main="")
  hist(thetas[(M0+1):niter],prob=TRUE,xlab=expression(theta),main="")
  th = seq(-6,6,length=1000)
  lines(th,dt(th,df=4),col=2,lwd=2)
}



# Target   = 0.1*N(0,1) + 0.9*N(10,1)
# Proposal = Uniform(theta(i-1)-sd,theta(i-1)+sd)
target = function(theta){
  0.1*dnorm(theta)+0.9*dnorm(theta,10,1)
}

M0     = 10000
M      = 10000
niter  = M0+M
thetas = rep(0,niter)
par(mfrow=c(3,3))
for (sd in c(1,5,10)){
  theta = -5
  for (iter in 1:niter){
    theta.star = runif(1,theta-sd,theta+sd)
    acceptance = min(1,target(theta.star)/target(theta))
    if (runif(1)<acceptance){
      theta = theta.star
    }
    thetas[iter] = theta
  }
  ts.plot(thetas,xlab="Iterations",ylab=expression(theta))
  abline(v=M0,col=2)
  thetas = thetas[(M0+1):niter]
  title(paste("SD=",sd,sep=""))
  acf(thetas,main="")
  hist(thetas,prob=TRUE,xlim=range(thetas,th),main="",xlab=expression(theta))
  th = seq(-5,15,length=1000)
  lines(th,target(th),col=2,lwd=2)
}