## Monday, July 18, 2011

### Minimalist MCMC, in R

Christian Robert recently posted a minimalist piece of R code implementing a Metropolis sampler for an arbitrary target function. Here's the code with a few edits:
 it <- 10^5 target = function(x) (x>-2)*(x<2)*dnorm(x=x, mean=0, sd=1) mcmc=rep(0,it) for ( t in 2:it ) {
   prop = mcmc[t-1] + runif(n = 1, min = -0.4, max = 0.4)    if ( runif(1) < target(prop) / target(mcmc[t-1]) ) {       mcmc[t] = prop    } else {       mcmc[t] = mcmc[t-1]    } } hist(mcmc) plot( x = 1:it, y = mcmc, pch = '*' ) 

This was just throw-away code used to answer a simple question about MCMC, but it shows that the code can be pretty minimal.  People writing their first MCMC samplers, even if they are familiar with programming often make a predictable series of mistakes and writing some example samplers which deal with these problems in this format might be a useful approach to guiding them through it.

Some of the easily avoided problems which trip people up:
1. Trying to work with densities, rather than log-densities.  It fails fairly quickly.
2. Failing to deal with boundaries of parameter space (the sampler above _does_ deal with them in the target function).
3. Ignoring obvious caching opportunities.
Those are probably enough goals for minimal MCMC code for now.

#### 1 comment:

1. Thanks: I always strive to code in a minimal(ist) way, much to the sorrow of my co-authors...