Jason Blevins’s Notebook
Chib and Greenberg (1995)

Citation

Chib, Siddhartha and Edward Greenberg (1995): “Understanding the Metropolis Hastings Algorithm,” American Statistical Journal, 49, 327–335.

@Article{chib95understanding,
  author       = {Siddhartha Chib and Edward Greenberg},
  title        = {Understanding the {Metropolis-Hastings} Algorithm},
  journal      = {American Statistical Association},
  year         = 1995,
  volume       = 49,
  pages        = {327--335}
}

Introduction

The Metropolis-Hastings algorithm is a Markov chain Monte Carlo method developed by Metropolis, Rosenbluth, Rosenbluth, Teller, and Teller (1953) and generalized by Hastings (1970). This algorithm is very general and gives rise to the Gibbs sampler as a special case. This article provides a self-contained intuitive development of the algorithm from first principles. The paper focuses on absolutely continuous target densities, but the algorithm also applies to discrete and mixed distributions.

Acceptance-Rejection Sampling

Classical sampling methods typically generate independent samples. One such method is the A-R method. The objective is to sample from an absolutely continuous target density $\pi (x)=f(x)/K$, where $x\in {ℝ}^d$, $f(x)$ is the unnormalized target density, and $K$ the potentially unknown normalizing constant. Suppose that we can sample from another density $h(x)$ and that there exists a constant $c$ such that $f(x)\le ch(x)$ for all $x$. To obtain a draw from $\pi$:

1. Draw a candidate $Z$ from $h$ and $u$ from $U(\mathrm{0,1})$.
2. If $u\le \frac{f(Z)}{\mathrm{ch}(Z)}$, return $Z$.
3. Otherwise, return to 1.

The expected number of iterations required to accept a draw is $c^{-1}$. To ensure efficiency, the optimal choice of $c$ is

Markov Chain Monte Carlo Simulation

Consider a Markov transition kernel $P(x,A)$ for $x\in {ℝ}^d$ and $A\in ℬ$, where $ℬ$ is the Borel $\sigma$-algebra on ${ℝ}^d$. Two major concerns in Markov chain theory are whether there exists an invariant distribution and whether iterating the transition kernel converges to the invariant distribution. The invariant distribution ${\pi }^{*}$ satisfies

where $\pi$ is the density of ${\pi }^{*}$ with respect to Lebesgue measure. Let

denote the $n$-th iteration of the transition kernel. Conditions are described under which $P^{(n)}(x,A)$ converges to ${\pi }^{*}$ as $n\to \mathrm{\infty }$. That is, the distribution of the draws generated by the iterations is approximately ${\pi }^{*}$.

MCMC methods look at the problem from the opposite perspective. The invariant distribution is known (the target distribution). The problem is how to generate an appropriate transition kernel with the aforementioned convergence property. Suppose that the transition kernel can be expressed as

where $p(x,x)=0$, ${\delta }_x(\mathrm{dy})=1$ if $x\in \mathrm{dy}$ and $0$ otherwise, and define $r(x)=1-{\int }_{{ℝ}^d}p(x,y)\mathrm{dy}$. The $r(x)$ term represents the probability of staying at $x$ which may be nonzero. If this kernel satisfies the reversibility condition

then $\pi$ is the invariant distribution of $P(x,\cdot )$. The Metropolis-Hastings algorithm generates a $p(x,y)$ with this property.

The Metropolis-Hastings Algorithm

Let $q(x,y)$ denote the candidate-generating density, where $\int q(x,y)\mathrm{dy}=1$. This is essentially the conditional density of $y$ given $x$. This density could potentially be used as the $p(x,y)$ term in the transition kernel, but it may not satisfy (1). If, for example, we have

then we can adjust $q$ by using a probability of move $\alpha (x,y)$. Transitions will be made using

for $x\ne y$.

The choice of $\alpha$ follows the following logic. If (2) holds, then moves from $x$ to $y$ are happening too often under $q$. We should thus choose $\alpha (y,x)=1$. But then, in order to satisfy (1), we must have

This implies that

Conversely, we can consider the case when the inequality in (1) were reversed to derive $\alpha (y,x)$. Thus, to summarize, in order to satisfy reversibility, we set

Hence, the desired transition kernel is

Thus, the Metropolis-Hastings algorithm is defined by the candidate-generating density $q(x,y)$. Note that $\alpha (x,y)$ does not require knowledge of the normalizing constant because it drops out of the ratio $\pi (y)/\pi (x)$. A special case arises when the candidate-generating density is symmetric: $q(x,y)=q(y,x)$ since the probability of move reduces to $\pi (y)/\pi (x)$. This special case forms the basis for optimization algorithms such as Simulated Annealing.

The algorithm proceeds as follows: Given an initial value $x^{(0)}$, for each $j=1,2,\dots ,N$

1. Draw $y$ from $q(x^{(j)},\cdot )$ and $u$ from $U(\mathrm{0,1})$
2. If $u\le \alpha (x^{(j)},y)$, set $x^{(j+1)}=y$.
3. Otherwise, set $x^{(j+1)}=x^{(j)}$.