Animating the Metropolis algorithm

September 8, 2013
By

(This article was first published on Ecology in silico, and kindly contributed to R-bloggers)

The Metropolis algorithm, and its generalization (Metropolis-Hastings algorithm) provide elegant methods for obtaining sequences of random samples from complex probability distributions. When I first read about modern MCMC methods, I had trouble visualizing the convergence of Markov chains in higher dimensional cases. So, I thought I might put together a visualization in a two-dimensional case.

I’ll use a simple example: estimating a population mean and standard deviation. We’ll define some population level parameters, collect some data, then use the Metropolis algorithm to simulate the joint posterior of the mean and standard deviation.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
# population level parameters
mu <- 7
sigma <- 3

# collect some data (e.g. a sample of heights)
n <- 50
x <- rnorm(n, mu, sigma)

# log-likelihood function
ll <- function(x, muhat, sigmahat){
  sum(dnorm(x, muhat, sigmahat, log=T))
}

# prior densities
pmu <- function(mu){
  dnorm(mu, 0, 100, log=T)
}

psigma <- function(sigma){
  dunif(sigma, 0, 10, log=T)
}

# 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 = .1){ # must be symmetric
  x + rnorm(1, 0, dist)
}

iter = 10000
chains <- 3
posterior <- array(dim = c(chains, 2, iter))
accepted <- array(dim=c(chains, iter - 1))

for (c in 1:chains){
  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
    theta_star <- c(jump(theta.post[1, t-1]), jump(theta.post[2, t-1]))
    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]
    }
  }
  posterior[c, , ] <- theta.post
}

Then, to visualize the evolution of the Markov chains, we can make plots of the chains in 2-parameter space, along with the posterior density at different iterations, joining these plots together using ImageMagick (in the terminal) to create an animated .gif:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
require(sm)
seq1 <- seq(1, 300, by=5)
seq2 <- seq(300, 500, by=10)
seq3 <- seq(500, iter, by=300)
sequence <- c(seq1, seq2, seq3)
length(sequence)

xlims <- c(4, 10)
ylims <- c(1, 6)

dir.create("metropolis_ex")
setwd("metropolis_ex")

png(file = "metrop%03d.png", width=700, height=350)
  for (i in sequence){
    par(mfrow=c(1, 2))
    plot(posterior[1, 1, 1:i], posterior[1, 2, 1:i],
         type="l", xlim=xlims, ylim=ylims, col="blue",
         xlab="mu", ylab="sigma", main="Markov chains")
    lines(posterior[2, 1, 1:i], posterior[2, 2, 1:i],
          col="purple")
    lines(posterior[3, 1, 1:i], posterior[3, 2, 1:i],
          col="red")
    text(x=7, y=1.2, paste("Iteration ", i), cex=1.5)
    sm.density(x=cbind(c(posterior[, 1, 1:i]), c(posterior[, 2, 1:i])),
               xlab="mu", ylab="sigma",
               zlab="", zlim=c(0, .7),
               xlim=xlims, ylim=ylims, col="white")
    title("Posterior density")
  }
dev.off()
system("convert -delay 15 *.png metrop.gif")
file.remove(list.files(pattern=".png"))

To leave a comment for the author, please follow the link and comment on his blog: Ecology in silico.

R-bloggers.com offers daily e-mail updates about R news and tutorials on topics such as: visualization (ggplot2, Boxplots, maps, animation), programming (RStudio, Sweave, LaTeX, SQL, Eclipse, git, hadoop, Web Scraping) statistics (regression, PCA, time series, trading) and more...



If you got this far, why not subscribe for updates from the site? Choose your flavor: e-mail, twitter, RSS, or facebook...

Comments are closed.