# R code for multivariate random-walk Metropolis sampling

I couldn’t find a simple R code for random-walk Metropolis sampling (the symmetric proposal version of Metropolis Hastings sampling) from a multivariate target distribution in arbitrary dimensions, so I wrote one. This is also my first R code.
It requires the package MASS to sample from the multivariate normal proposal distribution using the mvrnorm function. If you are using R on Windows, you can download the package zip for Windows from the link, and use Packages > Install package(s) from local zip files… from the GUI to install the package.
The reason I couldn’t write the code for a general Metropolis algorithm (i.e. for any arbitrary symmetric proposal distribution, not just normal) or a more general Metropolis-Hastings algorithm (with any arbitrary proposal distribution, symmetric or not) is that generating the proposal point would then require sampling from an arbitrary proposal distribution. This is only easy for a few standard distributions, but hard in general (which is the point of using such algorithms in the first place).

I. Function

The following is the function that does the Random Walk Metropolis-Hastings sampling when supplied with the required arguments. Notes about the arguments follow the code.

```rwmetro <- function(target,N,x,VCOV,burnin=0)
{
require(MASS)   #requires package MASS for normal sampling
samples <- x
for (i in 2:(burnin+N))
{
prop <- mvrnorm(n = 1, x, VCOV)
if (runif(1) < min(1, target(prop)/target(x)))
x <- prop
samples <- rbind(samples,x)
}
samples[(burnin+1):(N+burnin),]
}
```

II. Arguments

1. target function
The function defining the multivariate target distribution, written as a function of an n-vector, where n is the number of variables on which the distribution is defined. The different variables of your distribution must be written as `x[1], x[2]` etc.
Example
An example is the following, defining the function f(x,y)=exp(-5*abs(x^2+y^2-1)) in two dimensions:

```ring2D <- function(x)	# x is a vector
{
exp(-5*abs(x[1]^2+x[2]^2-1))
}
```
2. N integer
The final sample size (i.e., excluding the burn-in length).
Example
`100`
3. x numeric vector
The starting point (vector) of your Metropolis-Hastings chain. This vector needs to be the same length as the dimension of the target distribution. Your target function must be able to accept this vector as an argument.
Example
`c(-0.1,20.5)`
4. VCOV numeric matrix
The variance-covariance matrix for the multivariate normal that is used as the proposal distribution for random-walk Metropolis-Hastings. The length of this matrix also has to be the same as the dimension of the target distribution, i.e. the length of vectors acceptable by `target` and the length of `x`. You can vary the entries of this matrix and observe your results to see what works better for sampling your target distribution.
Example
The following line defines a matrix in two dimensions with .01 variance for each variable and no covariance between them.

```vcov2D <- .01*diag(2)
```
5. burnin (optional) integer
The ‘burn-in’ length for the chain. The number specified will be the number of initial samples chucked. If nothing is specified, it’s taken to be 0.
Example
`20`

III. Output
numeric matrix
The output is a matrix where each row is a sample from your target distribution, excluding the initial burn-in samples. The number of rows is thus the sample size, and the number of columns is equal to the dimension of the target distribution. You can use this matrix however you want, to save, visualize or analyze your sample.
Example
An example output of sampling a 2D distribution with sample size 5. Each row is an x,y sample.
``` [1,] 0.12961923 0.03708061 [2,] 0.10765053 -0.02798036 [3,] 0.01112930 -0.07255766 [4,] 0.06049625 -0.04546265 [5,] 0.1832164 -0.1244354 ```

IV. Usage Example
Let’s take the target distribution we used as an example above, f(x,y)=exp(-5*abs(x^2+y^2-1)). This looks like a ring of radius 1 rising from the x-y plane and centered at the origin. Here is a gnuplot surface plot of the distribution (because I found it frustratingly hard to figure out a level plot in R):

Let’s generate a sample of size 40,000 from this distribution with the starting point (0,0) and without any burn-in length, and with the variance-covariance matrix we defined before. This is done by calling the function with the correct arguments:

```ringsample2D <- rwmetro(ring2D,40000,c(0,0), vcov2D)
```

This assumes that you have already defined the target `ring2D` and the matrix `vcov2D` in the way explained, preceding this function call.
`ringsample` now contains a random sample from the distribution.
With the following one-line R code, I made a scatter plot of the sample from the `ringsample` matrix.

```plot(ringsample2D[,1], ringsample2D[,2], xlim=c(-1.5,1.5),ylim=c(-1.5,1.5), main='Metropolis-Hastings Sample',xlab='x', ylab='y', pch='.')
```

The following is the plot:

Putting all the code together in sequence, this is what the full code for defining the arguments, drawing the sample, and making the plot for this example looks like:

```# Define arguments
ring2D <- function(x)    # x is a vector
{
exp(-5*abs(x[1]^2+x[2]^2-1))
}

vcov2D <- .01*diag(2)

#Define the sampling function
rwmetro <- function(target,N,x,VCOV,burnin=0)
{
require(MASS)   #requires package MASS for normal sampling
samples <- x
for (i in 2:(burnin+N))
{
prop <- mvrnorm(n = 1, x, VCOV)
if (runif(1) < min(1, target(prop)/target(x)))
x <- prop
samples <- rbind(samples,x)
}
samples[(burnin+1):(N+burnin),]
}

# Call the sampling function with the arguments
ringsample2D <- rwmetro(ring2D,40000,c(0,0),vcov2D)

# Use the sample
plot(ringsample2D[,1], ringsample2D[,2], xlim=c(-1.5,1.5),ylim=c(-1.5,1.5), main='Metropolis-Hastings Sample',xlab='x', ylab='y', pch='.')
```

The following is another sample, this time for a 3-variable function, f(x,y,z)=exp(-5*abs(x^2+y^2+z^2-1)). This is just the 3D version of the previous function, and its density is mostly concentrated on a spherical shell of radius 1 and of some thickness, centered at (0,0,0). Notice how we never have to spare a thought about normalizing these functions to get the sample, which is one of the advantages of the Metropolis-Hastings algorithm.
The target function is now:

```ring3D <- function(x)	# x is a vector
{
exp(-5*abs(x[1]^2+x[2]^2+x[3]^2-1))
}
```

and the variance-covariance matrix is, say, .01 down the diagonal again:

```vcov3D <- .01*diag(3)
```

Let’s go for a sample of 20000 points after a burn-in of 20000, starting from (0,0,0). The function is then called as:

```ringsample3D <- rwmetro(ring3D,20000,c(0,0,0),vcov3D,20000)
```

`ringsample3D` is now a 3×1000 matrix, and I use the following to get a rotatable 3D plot of it (requires the rgl package):

```require(rgl)
plot3d(ringsample3D[,1], ringsample3D[,2], ringsample3D[,3], xlab='x',ylab='y', zlab='z',col='red', size=1)
```

I would have preferred a spinning 3D gif, but I’m too tired to work that out right now.
Share the knowledge.

## 8 thoughts on “R code for multivariate random-walk Metropolis sampling”

1. Thanks for your explanation but I have a doubt concerning your acceptance probability in the algorithm. Shouldn’t it be min(1, target(prop)/target(x)*q(prop)/q(x)) ? Where q is your proposal and is defined by your normal function in the line before (mvnorm…)

Like

1. The additional ratio you’re talking about is q(current | prop)/q(prop | current).
This is the random-walk Metropolis, where the proposal q is a normal. Since that’s a symmetric function, this ratio is just 1.

Like

2. gmas80

Nice post!!!
There is a small typo for the 2d plotting. It should be: plot(ringsample2D[,1], ringsample2D[,2], xlim=c(-1.5,1.5),ylim=c(-1.5,1.5), main=”Metropolis-Hastings Sample”,xlab=”x”, ylab=”y”, pch=’.’)

Like

1. Fixed; thanks.

Like

3. Hey please can someone help me with R codes for Metropolis, Metropolis-Hastings and Gibbs samplers for bivariate normal distribution. I am really a novice but need it for some comparative work which is rather urgent. Thanks a lot.

Like

4. To speed things up consider changing these two lines. Save you from computing the cholesky
samples <- matrix(nrow=burnin+N,ncol=length(x))
prop <- rnorm(2)%*%VCOV

Like

5. How can we get the mixing time? Thanks!

Like

6. Coşkun Kuş

Can you provide gnuplot surface plot command for f(x,y)=exp(-5*abs(x^2+y^2-1)).

Like