# This script contains the following functions to illustrate
# MCMC algorithms:
#
# metropolis.update - step-by-step illustration of Metropolis
# metropolis - long run of Metropolis
# gibbs.update - step-by-step illustration of Gibbs Sampler
# gibbs - long version of Gibbs
# plot.mcmc - Plot time series summaries of the chain
################################################################
# METROPOLIS ALGORITHM
################################################################
#
#Suppose we want to sample from p(theta|y) ~ N( mu , sig^2)
#Take proposal distribution as normal centred on current point:
# q(theta* |theta_t) ~ N(theta_t, phi)
# Step-by-step version to illustrate how the algorithm works
metropolis.update =function(start=1, n=10, sig.q=1, mu.p=0, sig.p=1) {
#set up vector theta[] for the samples
theta=rep(NA,n)
theta[1]=start
#count the number of accepted proposals
accept=0
#Set up plots
layout(matrix(c(0,1,1,0,0,1,1,0,2,2,2,2,2,2,2,2),4,4,byrow=T))
#Selection of points
plot (c(-5,5),c(0,0.5),type="n",main="",xlab="theta",ylab="density")
#Target density
curve(dnorm(x,mu.p,sqrt(sig.p)),add=T,lty=2,col="blue")
#Proposal density
curve(dnorm(x,theta[1],sqrt(sig.q)),add=T,col="red")
#current proposal
lines(c(theta[1],theta[1]),c(0,0.5),col="red")
#Time series plot so far:
plot(c(0,n),c(min(start,-3),max(start,3)),
xlab="iterations",ylab="theta",type="n")
lines(c(0,n),c(1.96,1.96),lty=2,col="blue")
lines(c(0,n),-c(1.96,1.96),lty=2,col="blue")
readLines(n=1)
#Main loop
for (t in 2:n) {
#pick a candidate value theta* from N(theta[t-1], sig.q)
theta.star = rnorm(1, theta[t-1], sqrt(sig.q))
#generate a random number between 0 and 1:
u = runif(1)
# calculate acceptance probability:
r = (dnorm(theta.star,mu.p,sig.p)) / (dnorm(theta[t-1],mu.p,sig.p) )
a=min(r,1)
#if u2) lines(density(theta[1:t]),col="green",lty=2)
#Time series plot so far:
plot(c(0,n),c(min(start,-3),max(start,3)),
xlab="iterations",ylab="theta",type="n")
lines(0:t,theta[1:(t+1)],col="red")
lines(c(0,n),c(1.96,1.96),lty=2,col="blue")
lines(c(0,n),-c(1.96,1.96),lty=2,col="blue")
#Wait for keypress
readLines(n=1)
#end loop
}
}
#Main version - run for a 'long time' to see convergence
metropolis=function (start=1, n=100, sig.q=1, mu.p=0, sig.p=1) {
#starting value theta_0 = start
#run for n iterations
#set up vector theta[] for the samples
theta=rep(NA,n)
theta[1]=start
#count the number of accepted proposals
accept=0
#Main loop
for (t in 2:n) {
#pick a candidate value theta* from N(theta[t-1], sig.q)
theta.star = rnorm(1, theta[t-1], sqrt(sig.q))
#generate a random number between 0 and 1:
u = runif(1)
# calculate acceptance probability:
r = (dnorm(theta.star,mu.p,sig.p)) / (dnorm(theta[t-1],mu.p,sig.p) )
a=min(r,1)
#if u 1.96:\t\t",mean(samples[i,j]>1.96)*100,"%\n")
}
if (p>1) cat("\nCorrelation:",cor(samples[i,1],samples[i,2]),"\n")
invisible(samples[i,])
}