A simple Metropolis sampler Example

Professor Darren Wilkinson of Stochastic Modelling 
  School of Mathematics & Statistics, Newcastle University 


A simple Metropolis sampler

Let's look at simulating from a normal with zero mean and unit variance using a 
Metropolis algorithm with uniform proposal distribution. 

Implementation in R
A function for the Metropolis sampler for this problem is given below. The chain 
is initialised at zero, and at each stage a U(-alpha,alpha) innovation is 
proposed. 

norm<-function (n, alpha) 
{
        vec <- vector("numeric", n)
        x <- 0
        vec[1] <- x
        for (i in 2:n) {
                innov <- runif(1, -alpha, alpha)
                can <- x + innov
                aprob <- min(1, dnorm(can)/dnorm(x))
                u <- runif(1)
                if (u < aprob) 
                        x <- can
                vec[i] <- x
        }
        vec
}

So, innov is a uniform random innovation and can is the candidate point. aprob 
is the acceptance probability. The decision on whether or not to accept is then 
carried out on the basis of whether or not a U(0,1) is less than the acceptance 
probability. We can test this as follows. 

normvec<-norm(10000,1)
par(mfrow=c(2,1))
plot(ts(normvec))
hist(normvec,30)
par(mfrow=c(1,1))

This shows a well mixing chain and reasonably normal distribution for the 
values. However, this was based on a choice of alpha of 1. Other values of alpha 
will not affect the stationary distribution, but will affect the rate of mixing 
of the chain. 

Implementation in C.

Here is a function for the simulation loop of this problem. This function has 
pointer arguments so that it can be called from R. 

void metrop(long *np, double *alphap, double *vec)
{
  long n,i;
  double alpha,x,can,u,aprob;
  gsl_rng *r = gsl_rng_alloc(gsl_rng_mt19937);
  n=*np;alpha=*alphap;
  x=0.0;
  vec[0]=x;
  for (i=1;i<n;i++)
    {
      can=x+gsl_ran_flat(r,-alpha,alpha);
      aprob=min(1.0, gsl_ran_ugaussian_pdf(can)/gsl_ran_ugaussian_pdf(x));
      u=gsl_ran_flat(r,0.0,1.0);
      if (u < aprob)
	x=can;
      vec[i]=x;
    }
}

This relies on a min function and some GSL functions. 
The full C code for this example is available here as metrop.c in a form 
intended for dynamic loading into R. 
Other implementations
Those who are interested may also wish to look at corresponding code for 
LISP-STAT: metrop.lsp,Python: metrop.py and Sather: metrop.sa. 
          
      d.j.wilkinson@ncl.ac.uk 
      http://www.staff.ncl.ac.uk/d.j.wilkinson/ 


Some notes: we probably should keep the actual acceptance rate.

Another example is to have target density of "scaled inverse chi square"

fden <- function(x, n, a) { x^(-n/2)*exp(-a/(2*x)) }

the normalizing constant is not needed.