## Thursday, September 12, 2013

### 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