Metropolis–Hastings algorithm

Section of Markov chains, Gibbs sampling, and Metropolis–Hastings, part of Bayesian statistics

Metropolis–Hastings algorithm

Type Section
Taught in Markov chains, Gibbs sampling, and Metropolis–Hastings
Part Bayesian statistics
Questions using it 15

Why another sampler is needed

Main article: Why another sampler is needed

Gibbs sampling is wonderfully direct when every full conditional is a distribution from which we can simulate. Often it is not. The logit-normal coin model in the logit-normal Metropolis–Hastings exercise, for example, has posterior kernel

f(\theta) =\theta^{x-1}(1-\theta)^{n-x-1} \exp\!\left[-\frac12\{\operatorname{logit}(\theta)-1\}^2\right], \qquad 0<\theta<1.

This is not a Beta kernel, nor another standard density from which we can draw directly. We can nevertheless evaluate f(\theta) at any admissible value. Writing

\pi(\theta)=\frac{f(\theta)}{C}, \qquad C=\int_0^1 f(t)\,\mathrm{d} t,

makes the computational problem precise: we want draws from \pi using pointwise evaluations of f, without first calculating the intractable normalising constant C.

The Metropolis–Hastings strategy is to construct a Markov chain whose stationary distribution is \pi. After warm-up, its dependent states can be used to estimate posterior quantities. The design problem is therefore not to sample from \pi in one step, but to devise a transition rule that preserves \pi.

Deriving the acceptance probability

Main article: Deriving the acceptance probability

Begin with any proposal density q(y\mid x) that is easy to simulate. It will not generally satisfy detailed balance for \pi:

\pi(x)q(y\mid x)\neq \pi(y)q(x\mid y).

We correct the imbalance by accepting a proposed move x\to y with probability a(x,y). For y\neq x, the resulting transition density is

q_{\mathrm{MH}}(y\mid x)=q(y\mid x)a(x,y).

Detailed balance will hold if

\pi(x)q(y\mid x)a(x,y) =\pi(y)q(x\mid y)a(y,x).

Consequently the paired acceptance probabilities must obey \frac{a(x,y)}{a(y,x)} =\frac{\pi(y)q(x\mid y)}{\pi(x)q(y\mid x)} =:r(x,y). Only their ratio is prescribed. Since rejection wastes a possible move, take the largest probabilities compatible with Equation (4) in Markov chains, Gibbs sampling, and Metropolis–Hastings and with 0\leq a\leq1. If r(x,y)\leq1, set a(x,y)=r(x,y) and a(y,x)=1; if r(x,y)>1, set a(x,y)=1 and a(y,x)=1/r(x,y). Both cases are a(x,y)=\min\!\left\{ 1, \frac{\pi(y)q(x\mid y)}{\pi(x)q(y\mid x)} \right\}.

Definition 1 (Metropolis–Hastings update) Let \pi(x) be a target density known only up to proportionality. Given X_t=x:

  1. propose Y=y from a proposal density q(y\mid x);

  2. calculate r(x,y) from Equation (4) in Markov chains, Gibbs sampling, and Metropolis–Hastings and put a(x,y)=\min\{1,r(x,y)\};

  3. draw U\sim\operatorname{Unif}(0,1) and set

X_{t+1}= \begin{cases} y,&U\leq a(x,y),\\ x,&U>a(x,y). \end{cases}

If \pi(x)=f(x)/C for unknown C, then

\frac{\pi(y)}{\pi(x)} =\frac{f(y)/C}{f(x)/C} =\frac{f(y)}{f(x)},

so the normalising constant cancels.

Article: Metropolis–Hastings update →

For a posterior, the evidence p(\text{data}) is precisely such a constant, so the ratio needs only likelihood times prior. In computation we normally evaluate the log ratio and accept when

\log U\leq \min\{0,\log r(x,y)\}.

This avoids numerical underflow when an unnormalised target is a product of many small terms.

This cancellation is why Metropolis methods became central in statistical physics. An Ising model on m binary sites has 2^m configurations, making its partition function impractical to sum, but ratios of unnormalised target weights are immediate. The original Metropolis algorithm used symmetric proposals; Hastings’ generalisation supplied the proposal-density correction in Equation (5) in Markov chains, Gibbs sampling, and Metropolis–Hastings.

The wartime history surrounding Los Alamos also includes Leó Szilárd.

For y\neq x, the Metropolis–Hastings transition density is

q_{\mathrm{MH}}(y\mid x)=q(y\mid x)a(x,y).

For y=x, the transition probability includes the total rejection probability:

q_{\mathrm{MH}}(x\mid x) =1-\int_{z\neq x}q(z\mid x)a(x,z)\,\mathrm{d} z.

Thus a rejected proposal is not thrown away in search of another proposal: the current value is recorded again. If instead we redraw until a proposal is accepted, the transition density becomes

\widetilde q(y\mid x) =\frac{q(y\mid x)a(x,y)}{A(x)}, \qquad A(x)=\int q(z\mid x)a(x,z)\,\mathrm{d} z.

The extra normaliser depends on the current state. This altered chain is in detailed balance with a density proportional to \pi(x)A(x), not generally with \pi(x), and therefore targets the wrong distribution.

Detailed-balance proof without skipped cases

Main article: Detailed-balance proof without skipped cases

Theorem 1 (Metropolis–Hastings invariance) The Metropolis–Hastings transition defined by Equation (5) in Markov chains, Gibbs sampling, and Metropolis–Hastings satisfies detailed balance with respect to \pi. Consequently, \pi is stationary.

Article: Metropolis–Hastings invariance →

Proof. Proof. Let

r(x,y)=\frac{\pi(y)q(x\mid y)} {\pi(x)q(y\mid x)}.

For y\neq x, \begin{align*} \pi(x)q_{\mathrm{MH}}(y\mid x) &=\pi(x)q(y\mid x)\min\{r(x,y),1\}. \end{align*} There are two cases.

If r(x,y)\leq1, then \begin{align*} \pi(x)q_{\mathrm{MH}}(y\mid x) &=\pi(x)q(y\mid x)r(x,y)\\ &=\pi(x)q(y\mid x) \frac{\pi(y)q(x\mid y)}{\pi(x)q(y\mid x)}\\ &=\pi(y)q(x\mid y). \end{align*} Since r(y,x)=1/r(x,y)\geq1, the reverse move is accepted with probability one, so

\pi(y)q_{\mathrm{MH}}(x\mid y) =\pi(y)q(x\mid y).

If r(x,y)>1, then the forward move is accepted with probability one:

\pi(x)q_{\mathrm{MH}}(y\mid x)=\pi(x)q(y\mid x).

Now r(y,x)<1, so \begin{align*} \pi(y)q_{\mathrm{MH}}(x\mid y) &=\pi(y)q(x\mid y)r(y,x)\\ &=\pi(y)q(x\mid y) \frac{\pi(x)q(y\mid x)}{\pi(y)q(x\mid y)}\\ &=\pi(x)q(y\mid x). \end{align*} Thus detailed balance holds in both cases. For x=y, both sides of the detailed-balance equation are identical. Therefore \pi is stationary. 

The theorem establishes invariance, not convergence from an arbitrary starting value and not efficiency. The ergodicity conditions stated earlier are still needed for convergence, while proposal design determines how much useful information the dependent output contains.

Used in

See also

Metropolis–Hastings update · Metropolis–Hastings invariance · Metropolis algorithm · Hastings correction · Markov chain · Stationarity and detailed balance · Detailed balance implies stationarity · Systematic-scan Gibbs sampler

Glossary and index