Markov Chain Monte Carlo
and Gibbs Sampling
Lecture Notes for EEB 581, version 26 April 2004
c
°B. Walsh 2004
A major limitation towards more widespread implementation of Bayesian ap-
proaches is that obtaining the posterior distribution often requires the integration
of high-dimensional functions. This can be computationally very difficult, but
several approaches short of direct integration have been proposed (reviewed by
Smith 1991, Evans and Swartz 1995, Tanner 1996). We focus here on Markov
Chain Monte Carlo (MCMC) methods, which attempt to simulate direct draws
from some complex distribution of interest. MCMC approaches are so-named be-
cause one uses the previous sample values to randomly generate the next sample
value, generating a Markov chain (as the transition probabilities between sample
values are only a function of the most recent sample value).
The realization in the early 1990’s (Gelfand and Smith 1990) that one particu-
lar MCMC method, the Gibbs sampler, is very widely applicable to a broad class
of Bayesian problems has sparked a major increase in the application of Bayesian
analysis, and this interest is likely to continue expanding for sometime to come.
MCMC methods have their roots in the Metropolis algorithm (Metropolis and
Ulam 1949, Metropolis et al. 1953), an attempt by physicists to compute com-
plex integrals by expressing them as expectations for some distribution and then
estimate this expectation by drawing samples from that distribution. The Gibbs
sampler (Geman and Geman 1984) has its origins in image processing. It is thus
somewhat ironic that the powerful machinery of MCMC methods had essentially
no impact on the field of statistics until rather recently. Excellent (and detailed)
treatments of MCMC methods are found in Tanner (1996) and Chapter two of
Draper (2000). Additional references are given in the particular sections below.
MONTE CARLO INTEGRATION
The original Monte Carlo approach was a method developed by physicists to use
random number generation to compute integrals. Suppose we wish to compute
a complex integral
Z
b
a
h(x) dx (1a)
If we can decompose h(x) into the production of a function f(x) and a probability
1
