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

A simple Metropolis-Hastings independence sampler

Let's look at simulating from a gamma distribution with arbitrary shape and scale parameters, using a Metropolis-Hastings independence sampling algorithm with normal proposal distribution with the same mean and variance as the desired gamma.

Implementation in R

A function for the Metropolis-Hastings sampler for this problem is given below. The chain is initialised at zero, and at each stage a N(a/b,a/(b*b)) candidate is proposed. 
gamm<-function (n, a, b) 
{
        mu <- a/b
        sig <- sqrt(a/(b * b))
        vec <- vector("numeric", n)
        x <- a/b
        vec[1] <- x
        for (i in 2:n) {
                can <- rnorm(1, mu, sig)
                aprob <- min(1, (dgamma(can, a, b)/dgamma(x, 
                        a, b))/(dnorm(can, mu, sig)/dnorm(x, 
                        mu, sig)))
                u <- runif(1)
                if (u < aprob) 
                        x <- can
                vec[i] <- x
        }
        vec
}
So, can is the candidate point and 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. 
vec<-gamm(10000,2.3,2.7)
par(mfrow=c(2,1))
plot(ts(vec))
hist(vec,30)
par(mfrow=c(1,1))
This shows a well mixing chain and reasonably gamma-like distribution for the values.

The full R source code for this example is available here as indep.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 gamm(long *np, double *ap, double *bp, double *vec)
{
  long n,i;
  double a,b,x,can,u,aprob,mu,sig;
  gsl_rng *r = gsl_rng_alloc(gsl_rng_mt19937);
  n=*np;a=*ap;b=*bp;
  mu=a/b;
  sig=sqrt(a/(b*b));
  x=mu;
  vec[0]=x;
  for (i=1;i<n;i++)
    {
      can=mu+gsl_ran_gaussian(r,sig);
      aprob=min(1.0, (gsl_ran_gamma_pdf(can,a,1/b)/gsl_ran_gamma_pdf(x,a,1/b))/
		(gsl_ran_gaussian_pdf(can-mu,sig)/
		 gsl_ran_gaussian_pdf(x-mu,sig)));
      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 indep.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: indep.lsp,Python: indep.py and Sather: indep.sa.

Darren J Wilkinson Book: Stochastic Modelling for Systems Biology  
darren.wilkinson@ncl.ac.uk
http://www.staff.ncl.ac.uk/d.j.wilkinson/