# 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.