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

By BioStatMatt

(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 

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 

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…

Source:: R News

Leave a Reply

Your email address will not be published. Required fields are marked *

Time limit is exhausted. Please reload CAPTCHA.