Prof Darren Wilkinson Professor 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.

The full R source code for this example is available here as metrop.r.

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.

Darren J Wilkinson

darren.wilkinson@ncl.ac.uk

http://www.staff.ncl.ac.uk/d.j.wilkinson/

	Professor of Stochastic Modelling
	School of Mathematics & Statistics
	Newcastle University


		`darren.wilkinson@ncl.ac.uk`
		`http://www.staff.ncl.ac.uk/d.j.wilkinson/`