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.