Quick illustration of Metropolis and Metropolis-in-Gibbs Sampling in R

June 3, 2017
By

(This article was first published on R – BioStatMatt, and kindly contributed to R-bloggers)

The code below gives a simple implementation of the Metropolis and Metropolis-in-Gibbs sampling algorithms, which are useful for sampling probability densities for which the normalizing constant is difficult to calculate, are irregular, or have high dimension (Metropolis-in-Gibbs).


## Metropolis sampling
## x    - current value of Markov chain (numeric vector)
## targ - target log density function
## prop - function with prototype function(x, ...) that generates 
##   a proposal value from a symmetric proposal distribution
library('mvtnorm')
prop_mvnorm <- function(x, ...)
  drop(rmvnorm(1, mean=x, ...))
metropolis <- function(x, targ, prop=prop_mvnorm, ...) {
  xnew <- prop(x)
  lrat <- targ(xnew, ...) - targ(x, ...)
  if(log(runif(1)) < lrat)
    x <- xnew
  return(x)
}

## Metropolis-in-Gibbs sampling
## x    - current value of Markov chain (numeric vector)
## targ - target log density function
## ...  - arguments passed to 'targ'
gibbs <- function(x, targ, ...) {
  for(i in 1:length(x)) {
    ## define full conditional
    targ1 <- function(x1, ...) {
      x[i] <- x1
      targ(x, ...)
    }
    ## sample using Metropolis algorithm
    x[i] <- metropolis(x[i], targ1, ...)
  }
  return(x)
}

The following code produces the figure below to illustrate the two methods using a 'dumbell' distribution (cf. R package 'ks' vignette).


### The code below illustrates the use of the functions above 

## target 'dumbell' density (c.f., R package 'ks' vignette)
library('ks')
mus <- rbind(c(-2,2), c(0,0), c(2,-2))
sigmas <- rbind(diag(2), matrix(c(0.8, -0.72, -0.72, 0.8), nrow=2), diag(2)) 
cwt <- 3/11
props <- c((1-cwt)/2, cwt, (1-cwt)/2)
targ_test <- function(x, ...)
  log(dmvnorm.mixt(x, mus=mus, Sigmas=sigmas, props=props))

## plot contours of target 'dumbell' density
set.seed(42)
par(mfrow=c(1,2))
plotmixt(mus=mus, Sigmas=sigmas, props=props,
         xlim=c(-4,4), ylim=c(-4,4),
         xlab=expression(x[1]),
         ylab=expression(x[2]),
         main="Metropolis-in-Gibbs")

## initialize and sample using Metropolis-in-Gibbs
xcur <- gibbs(c(0,0), targ_test, sigma=vcov_test)
for(j in 1:2000) {
  xcur <- gibbs(xcur, targ_test, sigma=vcov_test)
  points(xcur[1], xcur[2], pch=20, col='#00000055')
}

## plot contours of target 'dumbell' density
plotmixt(mus=mus, Sigmas=sigmas, props=props,
         xlim=c(-4,4), ylim=c(-4,4),
         xlab=expression(x[1]),
         ylab=expression(x[2]),
         main="Metropolis")

## initialize and sample using Metropolis
xcur <- metropolis(c(0,0), targ_test, sigma=vcov_test)
for(j in 1:2000) {
  xcur <- metropolis(xcur, targ_test, sigma=vcov_test)
  points(xcur[1], xcur[2], pch=20, col='#00000055')
}

The figure illustrates two contrasting properties of the two methods:

  1. Metropolis-in-Gibbs samples can get 'stuck' in certain regions of the support, especially when there are multiple modes and/or significant correlation among the random variables. This is not as much a problem for Metropolis sampling.
  2. Metropolis sampling can produce fewer unique samples due to the poor approximation of the proposal density to the target density. This occurs more often for high-dimensional target densities.

To leave a comment for the author, please follow the link and comment on their blog: R – BioStatMatt.

R-bloggers.com offers daily e-mail updates about R news and tutorials on topics such as: Data science, Big Data, R jobs, 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.

Sponsors

Never miss an update!
Subscribe to R-bloggers to receive
e-mails with the latest R posts.
(You will not see this message again.)

Click here to close (This popup will not appear again)