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):

Target Distribution in Gnuplot
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:
Metropolis-Hastings Sample, scatter 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)

Metropolis Sample

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

Advertisements

3D Tic-Tac-Toe

It was during a boring lecture at NCRA that Debsankha and I came up with the idea of this game. We were wondering how to pass the time, when he suggested we play a 4X4 tic-tac-toe. I had tried it by myself once, and it wasn’t too good, so I suggested a 3X3X3 version instead. We drew out the grid, and played it, and believe me it was fun.

Let me give a brief idea of how to play this.

The tic-tac-toe is essentially being played in a 3X3X3 grid. All the rules are same as in the normal tic-tac-toe, except that now you have many more lines which you can complete to score. You have all the edges, face and body diagonals etc. of the cube, along with the normal 2D tic-tac-toe lines on each face of the cube. This gives you a total of more than 40 lines (haven’t calculated it exactly yet), as opposed to only 8 for the normal tic-tac-toe.

But the trick is, you cannot play this on a 3D cube. You must draw three normal tic-tac-toe grids, corresponding to each ‘floor’ of the cube, like so:

 

3D3T-0

Take any to be the lowest floor, doesn’t really matter (as the rules of the game are reflection-symmetric). Play as in the normal tic-tac-toe, but don’t cross off the cells once you complete a line. This is because first, only lines that belong entirely in one floor can be crossed off in the old style. Second, and this is interesting, more than one of your lines may pass through a single cell. So instead, you must add to your count of lines. What’s more, before you add to your score, you must explain to your opponent how exactly you are completing a line. It may so happen that in one move you have completed more than one line, which you must also explain. But here’s the catch: if you fail to notice that you have completed a line, you will not report it, and hence miss that point. Here’s a 3D tic-tac-toe game in progress:

3D3T-1

I’ve pointed out the completed lines with the same colour. The circle with the red/blue gradient belongs to both the red and the blue line. The yellow and blue lines are verticals, for example, the green is a body diagonal, and the red is a normal line on the top floor. The black ones do not complete any line.

In a normal game there won’t be any of these colours, so you cannot keep track of which lines have already been drawn. That might make it difficult for you to diagnose whether or not your opponent has really completed a new line when they are claiming to. Fortunately, there’s no need to keep track of this. Just make sure that when a player claims to have completed a new line, that line must involve the last cell which they have filled. The reason for this is simple, figure it out.

We noticed that the centre cell of the middle floor is vital, because many lines pass through it. If a player places their symbol there, they’ll block many lines. This doesn’t matter so much in 2D tic-tac-toe, which isn’t too smart a game anyway, but here the benefit is increased a lot. This means the starting player has an unfair advantage, because they can block that cell. Moreover, since there are 27 cells in total, the starting player will have one more move than the opponent, which is also a little unfair.

To simultaneously take care of both the problems, we decided to make that cell ‘democratic’, meaning that we shall assume that each player has a symbol in that cell, and their lines can pass through it. This greatly increases the number of possible lines in the game and makes it more interesting. Also, now there are only 26 cells to fill, half of which is filled by each player.

Therefore, we shall put a little symbol in that cell that signifies that it has already been allocated. What better symbol to put there than the ‘trademark’?

3D3T-2

3T stands for Tic-Tac-Toe. It’s just like the centre white cubie of the Rubik’s cube or the Ace of Spades in a deck of cards. It’s special.

Note one more thing. Since now we have many lines to complete, we can even have three players at this game, and everyone will have some lines to claim for themselves. We actually did play with three players. But we needed a third symbol now. So just to increase the geek quotient, we used the initials of our surnames in Greek letters as our symbols. So we had δ, μ and ɣ. You can try it too, if you want onlookers to be awed.

We thought of extending this to even higher dimensions, but then we’ll be left only with a tensor-like index notation, and no good way to visualize the grid. I thought that’s what would make it interesting (what makes this one interesting is that you must visualize the 3D grid in your mind), but then I thought that it’s too complicated.

Overall, I think it’s a fun game, and a good exercise in mental visualization and spatial manipulations. It also requires brains, unlike the normal 2D 3T.