# World’s simplest Gibbs sampler illustration
# Posterior of (theta1, theta2) is BVN(y1, y2, 1, 1, rho)
# ‘Full conditionals’ are
# theta1 | theta2, y ~ Normal( y1 + rho*(theta2-y2), 1 – rho^2)
# theta2 | theta1, y ~ Normal( y2 + rho*(theta1-y1), 1 – rho^2)
rm(list=ls()); set.seed(20210311);
rho <- 0.8; y <- c(0, 0); # To code your own Gibbs sampler, write the update function(s), # and a BuildChain function, that calls the updates update <- function(other.theta, rho, y) { rnorm(1, mean = y + rho*(other.theta-y), sd = sqrt(1 - rho^2)) } # Only one update rule required for this simple example Build.chain <- function(inits, chain.length, rho, y) { chain <- matrix(rep(chain.length*2), chain.length, 2) chain[1,] <- inits other.theta <- inits[2] for(t in 2:chain.length) { other.theta <- update(other.theta, rho, y=y[2]) chain[t,1] <- other.theta other.theta <- update(other.theta, rho, y=y[1]) chain[t,2] <- other.theta } return(chain) } # Build 4 chains, starting at (+/-2.5, +/-2.5) chain1 <- Build.chain(c( 2.5, 2.5), 1000, rho=rho, y=y) chain2 <- Build.chain(c(-2.5, 2.5), 1000, rho=rho, y=y) chain3 <- Build.chain(c(-2.5,-2.5), 1000, rho=rho, y=y) chain4 <- Build.chain(c( 2.5,-2.5), 1000, rho=rho, y=y) # Illustrative plots op <- par(mfrow=c(1,3)) plot(chain1[1:10,], type="l", xlim=c(-4,4), ylim=c(-4,4), xlab="theta1", ylab="theta2", main="First 10 updates of 4 Gibbs samplers") lines(chain2[1:10,]); lines(chain3[1:10,]); lines(chain4[1:10,]); points(x=c(2.5,-2.5,-2.5,2.5), y=c(2.5,2.5,-2.5,-2.5), pch=15, cex=1.5, col="red") plot(chain1[1:200,], type="l", xlim=c(-4,4), ylim=c(-4,4), xlab="theta1", ylab="theta2", main="First 200 updates of 4 Gibbs samplers") lines(chain2[1:200,]); lines(chain3[1:200,]); lines(chain4[1:200,]); points(x=c(2.5,-2.5,-2.5,2.5), y=c(2.5,2.5,-2.5,-2.5), pch=15, cex=1.5, col="red") foo <- rbind(chain1[501:1000,], chain2[501:1000,], chain3[501:1000,], chain4[501:1000,]) plot(foo, xlab="theta1", ylab="theta2", xlim=c(-4,4), ylim=c(-4,4), pch=19, cex=.5, main="Final 500 draws of 4 Gibbs samplers") rm(foo) par(op)