<p><div><div><div><div data-show-faces="false" data-height="21" data-width="100" data-layout="button_count" data-send="true" data-href="https://www.r-bloggers.com/quick-illustration-of-metropolis-and-metropolis-in-gibbs-sampling-in-r/"></div></div><div></div></div><div></div></div>
<div>
(This article was first published on <strong><a href="http://biostatmatt.com/archives/2964"> R – BioStatMatt</a></strong>, and kindly contributed to <a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers)</a> </div>
<p>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).</p>
<pre><code>
## 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 &lt;- function(x, ...)
  drop(rmvnorm(1, mean=x, ...))
metropolis &lt;- function(x, targ, prop=prop_mvnorm, ...) {
  xnew &lt;- prop(x)
  lrat &lt;- targ(xnew, ...) - targ(x, ...)
  if(log(runif(1)) &lt; lrat)
    x &lt;- 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 &lt;- function(x, targ, ...) {
  for(i in 1:length(x)) {
    ## define full conditional
    targ1 &lt;- function(x1, ...) {
      x[i] &lt;- x1
      targ(x, ...)
    }
    ## sample using Metropolis algorithm
    x[i] &lt;- metropolis(x[i], targ1, ...)
  }
  return(x)
}
</code></pre>
<p>The following code produces the figure below to illustrate the two methods using a 'dumbell' distribution (cf. R package 'ks' vignette).</p>
<pre><code>
### The code below illustrates the use of the functions above 

## target 'dumbell' density (c.f., R package 'ks' vignette)
library('ks')
mus &lt;- rbind(c(-2,2), c(0,0), c(2,-2))
sigmas &lt;- rbind(diag(2), matrix(c(0.8, -0.72, -0.72, 0.8), nrow=2), diag(2)) 
cwt &lt;- 3/11
props &lt;- c((1-cwt)/2, cwt, (1-cwt)/2)
targ_test &lt;- 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=&quot;Metropolis-in-Gibbs&quot;)

## initialize and sample using Metropolis-in-Gibbs
xcur &lt;- gibbs(c(0,0), targ_test, sigma=vcov_test)
for(j in 1:2000) {
  xcur &lt;- 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=&quot;Metropolis&quot;)

## initialize and sample using Metropolis
xcur &lt;- metropolis(c(0,0), targ_test, sigma=vcov_test)
for(j in 1:2000) {
  xcur &lt;- metropolis(xcur, targ_test, sigma=vcov_test)
  points(xcur[1], xcur[2], pch=20, col='#00000055')
}
</code></pre>
<p><a rel="nofollow" href="https://i1.wp.com/biostatmatt.com/uploads/metropolis-in-gibbs.png" target="_blank"><img sizes="(max-width: 604px) 100vw, 604px" src="https://i1.wp.com/biostatmatt.com/uploads/metropolis-in-gibbs-625x349.png?w=450" alt="" class="wp-image-2965" data-recalc-dims="1"></a></p>
<p>The figure illustrates two contrasting properties of the two methods:</p>
<ol>
<li>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.</li>
<li>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.</li>
</ol>
	<p>
<div>
<div>To <strong>leave a comment</strong> for the author, please follow the link and comment on their blog: <strong><a href="http://biostatmatt.com/archives/2964"> R – BioStatMatt</a></strong>.</div>
<hr>
<a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers.com</a> offers <strong><a rel="nofollow" href="https://feedburner.google.com/fb/a/mailverify?uri=RBloggers">daily e-mail updates</a></strong> about <a rel="nofollow" href="https://www.r-project.org/" title="The R Project for Statistical Computing">R</a> news and <a rel="nofollow" href="https://www.r-bloggers.com/search/tutorial" title="R tutorials">tutorials</a> on topics such as: <a rel="nofollow" href="https://www.r-bloggers.com/search/data%20science" title="Data science">Data science</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/Big%20Data" title="Big Data">Big Data, <a rel="nofollow" href="https://www.r-users.com/" title="R jobs">R jobs</a>, visualization (<a rel="nofollow" href="https://www.r-bloggers.com/search/ggplot2" title="ggplot and ggplot2 tutorials">ggplot2</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/boxplot" title="Boxplots using lattice and ggplot2 tutorials">Boxplots</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/map" title="Maps and gis">maps</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/animation" title="Animation in R">animation</a>), programming (<a rel="nofollow" href="https://www.r-bloggers.com/search/RStudio" title="RStudio IDE for R">RStudio</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/sweave" title="Sweave and literate programming">Sweave</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/LaTeX" title="LaTeX in R">LaTeX</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/SQL" title="SQL and databases">SQL</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/eclipse" title="Eclipse IDE for R">Eclipse</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/git" title="git and github, Version Control System">git</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/hadoop" title="Large data in R using Hadoop">hadoop</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/Web+Scraping" title="Web Scraping of google, facebook, yahoo, twitter and more using R">Web Scraping</a>) statistics (<a rel="nofollow" href="https://www.r-bloggers.com/search/regression" title="Regressions and ANOVA analysis tutorials">regression</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/PCA" title="principal component analysis tutorial">PCA</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/time+series" title="Time series">time series</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/trading" title="finance trading">trading</a>) and more...

	<img width="1" alt="" src="http://feeds.feedburner.com/~r/RBloggers/~4/hnoMJ0v4Ufs" height="1"></a></div>