**Tags**

2D, 3D, code, gnuplot, Markov Chain, Metropolis Hastings, Monte Carlo, multi-dimensional, R, random walk, sampling, statistics

I couldn’t find a simple R code for random-walk 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. If you are using 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.

**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) { require(MASS) #requires package MASS samples <- x for (i in 2:N) { prop <- mvrnorm(n = 1, x, VCOV) if (runif(1) < min(1, target(prop)/target(x))) x <- prop samples <- cbind(samples,x) } samples }

**II. Arguments**

**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)) }

**N***integer*

Sample size

**Example**

`100`

**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)`

**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 <- matrix(data = c(.01,0,0,.01), nrow = 2, ncol = 2, byrow = TRUE)

**III. Output**

**samples***numeric matrix*

The output,`samples`

, is the final list of samples from your target distribution. It’s a matrix where each column is a sample. The number of rows is thus equal to the dimension of the target distribution, and the number of columns is the sample size. 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*column*is an x,y sample.

[1,] 0.12961923 0.10765053 0.01112930 0.06049625 0.1832164

[2,] 0.03708061 -0.02798036 -0.07255766 -0.04546265 -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 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='.')

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 <- matrix(data = c(.01,0,0,.01), nrow=2, ncol=2, byrow = TRUE) #Define the sampling function rwmetro <- function(target,N,x,VCOV) { require(MASS) #requires package MASS samples <- x for (i in 2:N) { prop <- mvrnorm(n = 1, x, VCOV) if (runif(1) < min(1, target(prop)/target(x))) x <- prop samples <- cbind(samples,x) } samples } # 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 <- matrix(data = c(.01,0,0,0,.01,0,0,0,.01), nrow=3, ncol=3, byrow = TRUE)

Let’s go for a sample of 40000 points again, starting from (0,0,0). The function is then called as:

ringsample3D <- rwmetro(ring3D,40000,c(0,0,0),vcov3D)

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