rm(list=ls()); set.seed(20210309);
# We reproduce the simple example describe on page 278 of
# Gelman et all
Build.chain <- function(inits, chain.length) { chain <- matrix(rep(chain.length*2), chain.length, 2) chain[1,] <- inits theta <- inits for(t in 2:chain.length) { theta.star <- rnorm(2, mean=theta, sd=1/5) log.r <- sum(dnorm(theta.star, 0, 1, log=T)) - sum(dnorm(theta, 0, 1, log=T)) if(log(runif(1)) < log.r){ theta <- theta.star } chain[t,] <- theta } return(chain) } # Construct 5 Markov chains from 5 different starting points chain1 <- Build.chain(inits=c( 3, 3), chain.length=1000) chain2 <- Build.chain(inits=c( 3,-3), chain.length=1000) chain3 <- Build.chain(inits=c(-3,-3), chain.length=1000) chain4 <- Build.chain(inits=c(-3, 3), chain.length=1000) chain5 <- Build.chain(inits=c( 0, 0), chain.length=1000) # Plot the results op <- par(mfrow=c(1,3)) # Plot (a): First 50 iterations sel <- 1:50 plot(jitter(chain5[sel,]), xlim=c(-4,4), ylim=c(-4,4), xlab="", ylab="", type="l", main="50 iterations"); lines(jitter(chain1[sel,])) ; lines(jitter(chain2[sel,])) ; lines(jitter(chain3[sel,])) ; lines(jitter(chain4[sel,])) ; points(x=c(3,3,-3,-3,0), y=c(3,-3,-3,3,0), col="red", pch=15, cex=1.5) rm(sel) # Plot (b): First 500 sel <- 1:500 plot(jitter(chain5[sel,]), xlim=c(-4,4), ylim=c(-4,4), xlab="", ylab="", type="l", main="500 iterations"); lines(jitter(chain1[sel,])) ; lines(jitter(chain2[sel,])) ; lines(jitter(chain3[sel,])) ; lines(jitter(chain4[sel,])) ; points(x=c(3,3,-3,-3,0), y=c(3,-3,-3,3,0), col="red", pch=15, cex=1.5) rm(sel) # Plot (c): Scatterplot of last 500, all combined sel <- 501:1000 foo <- rbind(chain1[sel,], chain2[sel,], chain3[sel,], chain4[sel,], chain5[sel,]) plot(jitter(foo), xlim=c(-4,4), ylim=c(-4,4), xlab="", ylab="", pch=19, cex=.5, main="5 chains combined") rm(foo); rm(sel); par(op)