Markov Chain Monte Carlo Simulation Using the Metropolis Algorithm

Manipulate[
SeedRandom[sr];
n=1000;
x1=-4;
y1=9;
xy=Table[{0,0},{i,1,n}];
xy[[1,1]]=x1;
xy[[1,2]]=y1;
na=0;
xypts=Table[xm=xy[[All,1]][[i-1]];ym=xy[[All,2]][[i-1]];ndistx=NormalDistribution[xm,proposalSigma];ndisty=NormalDistribution[ym,proposalSigma];
ux=RandomReal[ndistx];
uy=RandomReal[ndisty];
p2=pdf/.{x->ux,y->uy};p1=pdf/.{x->xm,y->ym};
r=p2/p1;u2=RandomReal[];
If[u2≤r,xy[[i]]={ux,uy};na=na+1,xy[[i]]={xm,ym}];
xy[[i]],
{i,2,n}];plgraphics=Graphics[Text[Style["proposal σ = "<>ToString[proposalSigma]<>"\nacceptance rate = "<>ToString[100.*na/n]<>"%",16,Bold],{-3,9},{-1,0},{1,0}]];
points=ListPlot[xypts,PlotRange->{{-5,10},{-5,10}},Frame->True,FrameLabel->{Style["x",Italic],Style["y",Italic]}];xpoints=ListPlot[xy[[All,1]],PlotRange->All,Frame->True,PlotLabel->Row[{Style["x",Italic]," trace"}],FrameLabel->{"iteration",Style["x",Italic]},ImagePadding->{{35,10},{35,25}}];ypoints=ListPlot[xy[[All,2]],PlotRange->All,Frame->True,PlotLabel->Row[{Style["y",Italic]," trace"}],FrameLabel->{"iteration",Style["y",Italic]},ImagePadding->{{35,10},{35,25}}];
GraphicsGrid[{{plgraphics,Show[plcont,points]},{xpoints,ypoints}}],{{proposalSigma,.2,"proposal σ"},.2,10.2,1.0},
{{sr,123,"random seed"},1,123456,1},TrackedSymbols:>Manipulate,
Initialization:>
{
pdf=0.087 (E^((1/2) (-(x-4) (0.595(x-4)-0.238 (y-3))-(-0.238(x-4)+0.595 (y-3))  (y-3))))+((E^((1/2) (-(x^2)-(y^2))))/(2 π));plcont=ContourPlot[pdf,{x,-4,10},{y,-4,10},PlotRange->{{-5,10},{-5,10},{0,0.25}},Contours->15,ContourShading->False,FrameLabel->{Style["x",Italic],Style["y",Italic]}];
}]

This Demonstration allows a simple exploration of the Metropolis algorithm sampling of a two-dimensional target probability distribution as a function of the width (sigma) of a circularly symmetric Gaussian proposal distribution. If the sigma is too small almost all proposals are accepted, the samples are highly correlated, and the burn-in period is long. If  the proposal sigma is too large, hardly any proposals are accepted, and the sampling sticks. When a proposal is not accepted the new sample is just a repeat of the previous sample. Somewhere in between, the samples are minimally correlated and efficiently concentrate in regions where the target probability distribution is significant.

The Demonstration produces three graphs. In the first the dots show the individual Metropolis samples in (x,y). The true target probability distribution is shown by the blue contours. The proposal sigma is given on the left along with the proposal acceptance rate expressed as a percentage. The second and third graphs show the individual samples in the x and y parameters. For small values of the proposal sigma, the burn-in period is sometimes longer than the simulation.

Monte Carlo methods provide approximate solutions to a great variety of problems in science and economics by performing statistical sampling experiments on a computer. The method applies to problems with no probabilistic content as well as to those with inherent probabilistic structure. Among all Monte Carlo methods, Markov chain Monte Carlo (MCMC) provides the greatest scope for dealing with very complicated systems. MCMC was first introduced in the early 1950s by statistical physicists (N. Metropolis, A. Rosenbluth, M. Rosenbluth, A. Teller, and E. Teller) as a method for the simulation of simple fluids. In the 1990s, the method began to play an important role in the arena of bioinformatics (the science of developing computer databases and algorithms to facilitate and expedite biological research, particularly in genomics).  The current renaissance in Bayesian statistics stems in large part on the use of MCMC methods to evaluate integrals in many dimensions. The earliest and still widely used MCMC method is called the Metropolis algorithm. In this Demonstration the Metropolis algorithm is used draw samples from a two-dimensional probability distribution in parameters, which consists of two elliptical Gaussians. Of course, the real value of the algorithm is in dealing with much higher-dimensional problems. The amazing property of the Metropolis algorithm is that after an initial burn-in period (which is discarded) the algorithm generates an equilibrium distribution of samples with a density proportional to the underlying probability distribution. It concentrates samples to regions with significant probability. These samples can then be used to obtain approximate integrals that are needed in many statistical problems.

Contributed by: Philip Gregory (Physics and Astronomy, University of British Columbia)