Visualization of Metropolis with rejected proposals

In response to a previous post, Tinu Schneider has embellished the Metropolis visualization of Maxwell Joseph so it includes explicit plotting of rejected proposals (shown as grey segments). Here are the animated gif and R code. Many thanks to Tinu for sending this! (The gif, below, may take a few moments to load.)

 
 
#
# Animating the Metropolis Algorithm
# -------------------------------------------
# source:
# http://mbjoseph.github.io/blog/2013/09/08/metropolis/
 
# Prepare
----------------------------------------
require(sm) # needs 'rpanel'
dir.create("metropolis_ex")
 
 
# Setup
# --------------------------------------
# population level parameters
mu    <- 6
sigma <- 3
 
# data-values
n <- 50 
 
# metropolis
iter   <- 10000
chains <-     3
sigmaJump <- 0.2
 
# plotting
seq1 <- seq(  1, 200,  by=  5)
seq2 <- seq(200, 450,  by= 10)
seq3 <- seq(450, 750,  by= 50)
seq4 <- seq(750, iter, by=150)
 
xlims <- c(4, 8)
ylims <- c(1, 5)
zlims <- c(0, 0.7)
 
widthPNG  <- 1200
heightPNG <-  750
 
 
 
# Helpers
# --------------------------------------
# collect some data (e.g. a sample of heights)
gendata <- function(n, mu, sigma) {
    rnorm(n, mu, sigma)
}
 
# log-likelihood function
ll <- function(x, muhat, sigmahat){
    sum(dnorm(x, muhat, sigmahat, log=TRUE))
}
 
# prior densities
pmu <- function(mu){
    dnorm(mu, 0, 100, log=TRUE)
}
 
psigma <- function(sigma){
    dunif(sigma, 0, 10, log=TRUE)
}
 
# posterior density function (log scale)
post <- function(x, mu, sigma){
    ll(x, mu, sigma) + pmu(mu) + psigma(sigma)
}
 
geninits <- function(){
    list(mu    = runif(1, 4, 10),
         sigma = runif(1, 2,  6))
}
 
jump <- function(x, dist = 0.1){ # must be symmetric
    x + rnorm(1, 0, dist)
}
 
 
# compute
# ----------------------------------------
 
 
# init 
proposed  <- array(dim = c(chains, 2, iter))
posterior <- array(dim = c(chains, 2, iter))
accepted  <- array(dim = c(chains, iter-1))
 
# the data
x <- gendata(n, mu, sigma)
 
for (c in 1:chains){
    props      <- array(dim=c(2, iter))
    theta.post <- array(dim=c(2, iter))
    inits <- geninits()
    theta.post[1, 1] <- inits$mu
    theta.post[2, 1] <- inits$sigma
    for (t in 2:iter){
        # theta_star = proposed next values for parameters
        jumpMu    <- jump(theta.post[1, t-1], sigmaJump)
        jumpSigma <- jump(theta.post[2, t-1], sigmaJump)
        theta_star <- c(jumpMu, jumpSigma)
        pstar <- post(x, mu = theta_star[1],      sigma = theta_star[2])
        pprev <- post(x, mu = theta.post[1, t-1], sigma = theta.post[2, t-1])
        lr <- pstar - pprev
        r  <- exp(lr)
 
        # theta_star is accepted if posterior density is higher w/ theta_star
        # if posterior density is not higher, it is accepted with probability r
        # else theta does not change from time t-1 to t
        accept <- rbinom(1, 1, prob = min(r, 1))
        accepted[c, t - 1] <- accept
        if (accept == 1){
            theta.post[, t] <- theta_star
        } else {
            theta.post[, t] <- theta.post[, t-1]
        }
        props[, t] <- theta_star
    }
    proposed[c, , ]  <- props
    posterior[c, , ] <- theta.post
}
 
 
# Plot
# ----------------------------------------------
 
oldPath <- getwd()
setwd("metropolis_ex")
 
sequence <- c(seq1, seq2, seq3, seq4)
cols <- c("blue", "purple", "red")
 
png(file = "metrop%03d.png", width=widthPNG, height=heightPNG)
    for (i in sequence){
        par(font.main=1, cex.main=1.2)
        layout(matrix(1:2, nrow=1), widths=c(1.0, 1.5))
  # prepare empty plot
        plot(0, type="n", 
             xlim=xlims, ylim=ylims,
             xlab="mu",  ylab="sigma", main="Markov chains")
        for (j in 1:chains) { # first all the gray lines
            lines(proposed[j, 1, 1:i],  proposed[j, 2, 1:i],  col="gray")
        }
        for (j in 1:chains) { # now the accepteds
            lines(posterior[j, 1, 1:i], posterior[j, 2, 1:i], col=cols[j])
        }
        text(x=6, y=1.0, paste("Iteration ", i), cex=1.2)
 
        sm.density(x=cbind(c(posterior[, 1, 1:i]), c(posterior[, 2, 1:i])),
                   xlab="mu",  ylab="sigma", zlab="", 
                   xlim=xlims, ylim=ylims,   zlim=zlims, 
                   col="white")
        title(main="Posterior density")
    }
dev.off()
 
 
# combine .png's into one .gif
# system('convert2gif.bat') # had to write a .bat-file
# system('convert -delay 20 *.png "metropolis.gif"') # doesent work
 
#file.remove(list.files(pattern=".png"))
 
setwd(oldPath)
 
#

Created by Pretty R at inside-R.org