# 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,]) }