Metropolis-Hastings

Metropolis-Hastings

Metropolis-Hastings (MH) is an MCMC algorithm that produces a sequence ${\theta }^{(0)},{\theta }^{(1)},\dots$ whose long-run distribution is a target density $\pi$ known only up to a normalizing constant. It works by converting an arbitrary proposal transition $Q(x,y)$ into a Markov chain that satisfies detailed balance with respect to $\pi$, using a single accept/reject step.

Given the current state ${\theta }^{(t)}=x$, one iteration proceeds as follows:

1. Propose $y\sim Q(x,\cdot )$.
2. Compute the acceptance probability
   $\alpha (x,y)=\min \left(1,\frac{\pi (y)Q(y,x)}{\pi (x)Q(x,y)}\right).$
3. Draw $u\sim \mathrm{Unif}(0,1)$. If $u<\alpha (x,y)$, accept: set ${\theta }^{(t+1)}=y$. Otherwise reject: set ${\theta }^{(t+1)}=x$.

The resulting chain $\{{\theta }^{(t)}\}$ satisfies detailed balance with respect to $\pi$:
$\pi (x)P(x,y)=\pi (y)P(y,x),$
where $P$ is the transition kernel of the MH chain. Under irreducibility, $\pi$ is therefore the unique stationary distribution, and by the ergodic theorem
$\frac{1}{N}{\sum }_{t=1}^{N}g({\theta }^{(t)})\stackrel{\text{a.s.}}{\to }\int g(\theta )\pi (\theta )d\theta$
for any integrable $g$. Only the ratio $\pi (y)/\pi (x)$ appears in $\alpha$, so any unknown normalizing constant of $\pi$ cancels. This is the reason MH is the workhorse of Bayesian computation: the unnormalized posterior $\tilde{\pi }(\theta )=\text{prior}(\theta )\text{likelihood}(\theta \mid \text{data})$ is sufficient to run the algorithm, and the intractable marginal likelihood is never needed.

When the proposal is symmetric — $Q(x,y)=Q(y,x)$, as when $y=x+Z$ for $Z$ symmetric about $0$ — the Hastings correction $Q(y,x)/Q(x,y)$ equals one and the acceptance probability collapses to
$\alpha (x,y)=\min \left(1,\frac{\pi (y)}{\pi (x)}\right).$
This special case is the Metropolis algorithm (proposed first historically; Hastings’ extension handles asymmetric $Q$). The further specialization in which $y=x+Z$ for a random-walk increment $Z$ is the random walk Metropolis chain, typically with $Z\sim N(0,\Sigma )$ or $Z\sim \mathrm{Unif}(-a,a{)}^d$.

Intuition

Standard rejection sampling requires a proposal density that uniformly dominates the target by a known constant — infeasible when $\pi$ is known only up to a normalizing constant, as happens for essentially every nontrivial Bayesian posterior. MH sidesteps the dominance requirement by constructing a Markov chain whose stationary distribution is $\pi$, rather than drawing independent samples outright. The chain’s samples are correlated, but long-run averages still converge to integrals against $\pi$ by the ergodic theorem.

The algorithm’s central trick is that an arbitrary proposal $Q$ will almost certainly produce a chain that fails detailed balance with respect to $\pi$. Adding a single accept/reject step with probability $\alpha (x,y)$ repairs the imbalance exactly. The ratio $\pi (y)Q(y,x)/\pi (x)Q(x,y)$ measures how much $Q$ overproposes the move $x\to y$ relative to what $\pi$ would like; the accept step discards the excess. The central tension is between proposals large enough to explore the state space — but so large they are rejected in the tails — and proposals small enough to be accepted but so small the chain crawls.

Properties

(Detailed Balance) The MH transition kernel satisfies $\pi (x)P(x,y)=\pi (y)P(y,x)$ by construction.

For $y\ne x$, the MH transition density factors as $P(x,y)=Q(x,y)\alpha (x,y)$. Multiplying by $\pi (x)$ and substituting the acceptance formula,
$\pi (x)P(x,y)=\pi (x)Q(x,y)\min \left(1,\frac{\pi (y)Q(y,x)}{\pi (x)Q(x,y)}\right)=\min (\pi (x)Q(x,y),\pi (y)Q(y,x)),$
and the right-hand side is symmetric under the swap $x\leftrightarrow y$. Detailed balance implies $\pi$ is stationary for $P$; irreducibility (plus aperiodicity in discrete state spaces) upgrades stationarity to convergence via the ergodic theorem. The full derivation is given in the Derivation section below.

(Only Ratios Needed) The acceptance probability depends on $\pi$ only through $\pi (y)/\pi (x)$, so unknown normalizing constants cancel.

If $\pi (\theta )=\tilde{\pi }(\theta )/Z$ for any constant $Z$, then $\pi (y)/\pi (x)=\tilde{\pi }(y)/\tilde{\pi }(x)$. This makes MH the standard tool for Bayesian inference: the unnormalized posterior $\tilde{\pi }(\theta )=\text{prior}(\theta )\text{likelihood}(\theta \mid \text{data})$ is enough to run the algorithm, and the marginal likelihood
$Z=\int \text{prior}(\theta )\text{likelihood}(\theta \mid \text{data})d\theta$
— the quantity that makes direct posterior sampling intractable — never has to be computed.

(Symmetric Proposals) If $Q(x,y)=Q(y,x)$, the Hastings correction drops out and $\alpha (x,y)=\min (1,\pi (y)/\pi (x))$.

The simplest instance is the random walk Metropolis chain with $y=x+Z$ for $Z$ symmetric about $0$ (e.g., $Z\sim N(0,\Sigma )$ or $Z\sim \mathrm{Unif}(-a,a{)}^d$). Under a symmetric $Q$, the accept step becomes: always accept moves to higher density; accept moves to lower density with probability equal to the density ratio. This preserves frequent exploration of high-density regions while still allowing downhill moves into the tails.

(Independence Sampler) If $Q(x,y)=q(y)$ does not depend on $x$, MH becomes an independence sampler with $\alpha (x,y)=\min (1,\frac{\pi (y)q(x)}{\pi (x)q(y)})$.

The proposal draws candidates from a fixed distribution $q$ independent of the current state. The Hastings correction $q(x)/q(y)$ is essential here — omitting it would target $q$, not $\pi$. Independence samplers mix well when $q$ is a close approximation to $\pi$ (a common choice is a Laplace approximation of the posterior, i.e., a Gaussian centered at the MAP with covariance $(-H\ell (\hat{\beta }){)}^{-1}$). They mix badly when $q$ has lighter tails than $\pi$, because the chain gets trapped in tail regions where $\pi /q$ is large and any move away is rejected.

(Proposal Scale Trade-off) The proposal variance controls an exploration-versus-acceptance trade-off; either extreme mixes slowly.

For random walk Metropolis with $Z\sim N(0,{\sigma }^{2}I)$: if $\sigma$ is too small, nearly every proposal is accepted but each move is tiny, producing strong autocorrelation and slow exploration of $\pi$. If $\sigma$ is too large, most proposals land in low-density regions and are rejected, producing a chain that rarely moves. The acceptance rate (fraction of proposals accepted) is the standard diagnostic: theory for high-dimensional Gaussian targets suggests tuning toward roughly $0.234$, and $0.2$–$0.4$ is a practical working range. When coordinates of $\theta$ live on very different scales, an isotropic ${\sigma }^{2}I$ proposal is the wrong shape — a diagonal or full covariance matched to the posterior (often obtained from a Laplace approximation of the posterior at the MAP) improves mixing dramatically.

(Burn-in and Initialization) Early iterates depend on ${\theta }^{(0)}$; they are discarded, and initializing near a posterior mode shortens the required burn-in.

The ergodic theorem guarantees $\frac{1}{N}\sum g({\theta }^{(t)})\to \int gd\pi$ as $N\to \infty$, but for finite $N$ the approximation depends on how quickly the chain forgets its starting point. Standard practice is to discard the first $B$ iterates — the burn-in — and use only ${\theta }^{(B+1)},\dots ,{\theta }^{(N)}$ for averages. Initializing at the MLE or MAP (rather than at an arbitrary point) cuts burn-in sharply, because the chain starts in a high-density region of $\pi$. Trace plots of each coordinate against iteration number are the standard visual diagnostic: a well-mixed chain looks like stationary noise around the posterior mean; a poorly-mixed chain drifts for many iterations or stays stuck at constant values.

(Any Proposal Works) Correctness of MH depends only on $\pi$, $Q$, and the acceptance formula; the choice of $Q$ determines efficiency, not whether the chain targets $\pi$.

Irreducibility requires only that $Q$ can reach every region of positive $\pi$-mass, however slowly. Different choices of $Q$ give qualitatively different algorithms — random walk Metropolis, independence samplers, Langevin dynamics, Hamiltonian Monte Carlo — but the underlying detailed-balance construction is identical. Choosing $Q$ is an engineering problem about mixing speed, not correctness. The acceptance formula is the unique maximally-accepting rule that preserves detailed balance: any smaller $\alpha$ that satisfies the balance equation could be enlarged toward $\min (1,\cdot )$ without breaking stationarity.

(Relation to Gibbs) Gibbs sampling is the special case of MH in which proposals are drawn from full conditionals and always accepted.

If $Q$ proposes ${\theta }^{\prime }$ by redrawing a single coordinate from its full conditional $\pi ({\theta }_j\mid {\theta }_{-j})$ and leaving the rest fixed, algebraic cancellation gives $\alpha (\theta ,{\theta }^{\prime })\equiv 1$, so every proposal is accepted. Gibbs thus avoids the acceptance overhead of MH but requires the full conditionals to be available in closed form. MH trades acceptance for the ability to handle arbitrary targets — the two algorithms sit at opposite ends of the same MCMC family, and hybrid “Metropolis-within-Gibbs” schemes handle cases where some conditionals are tractable and others are not.

Examples

Sampling from a non-standard 1D target

Suppose we want samples from $\pi (\theta )\propto \exp (-{\theta }^4/4)$ on $\mathbb{R}$, a target with no closed-form sampler. Use the random walk proposal $y=x+Z$ with $Z\sim N(0,{\sigma }^2)$.

Because $Q$ is symmetric, the acceptance probability collapses to
$\alpha (x,y)=\min \left(1,\frac{\exp (-y^4/4)}{\exp (-x^4/4)}\right)=\min (1,\exp (-(y^4-x^4)/4\left)\right).$

One iteration:

1. Draw $Z\sim N(0,{\sigma }^2)$ and set $y=x+Z$.
2. Compute $\log \alpha =\min (0,-(y^4-x^4)/4)$ (working in log-space avoids overflow when the exponent is large).
3. Draw $u\sim \mathrm{Unif}(0,1)$. If $\log u<\log \alpha$, accept (${\theta }^{(t+1)}=y$); otherwise reject (${\theta }^{(t+1)}=x$).

Run for $N=10,000$ iterations, discard the first $1,000$ as burn-in, and use the remaining $9,000$ samples to approximate integrals against $\pi$. With $\sigma \approx 1$ the acceptance rate lands near $50\%$ and the trace plot shows stationary fluctuation around the mode at $0$. With $\sigma =10$ almost every proposal is rejected and the trace plot is flat; with $\sigma =0.01$ almost every proposal is accepted but the chain drifts too slowly to cover the support in $N$ iterations.

Bayesian Poisson regression with a flat prior

The Poisson regression posterior with a flat prior on $\beta \in {\mathbb{R}}^{p+1}$ and observations $(y_i,x_i{)}_{i=1}^n$ is
$\pi (\beta \mid \text{data})\propto \exp \left({\sum }_{i=1}^n(y_i{x}_i^{\top }\beta -e^{x_i^{\top }\beta })\right),$
which is not a standard distribution. Propose updates from $Q(\beta ,\cdot )=N(\beta ,{\Sigma }_{\text{prop}})$ with ${\Sigma }_{\text{prop}}$ taken from the Laplace approximation at the MLE: ${\Sigma }_{\text{prop}}=c\cdot (-H\ell (\hat{\beta }){)}^{-1}$, where $\hat{\beta }$ is the MLE (from Newton’s method), $H\ell$ is the Hessian of the log-likelihood, and $c>0$ is a scaling factor tuned to hit the $0.2$–$0.4$ acceptance range.

Because $Q$ is symmetric (Gaussian centered at the current state), the Hastings correction drops out and
$\log \alpha (\beta ,{\beta }^{\prime })=\min \left(0,{\sum }_{i=1}^n(y_i{x}_i^{\top }({\beta }^{\prime }-\beta )-e^{x_i^{\top }{\beta }^{\prime }}+e^{x_i^{\top }\beta })\right).$

Initialize at $\hat{\beta }$ to skip most of the burn-in — initializing at $0$ instead puts the chain in a region where the Poisson mean $e^{x_i^{\top }\beta }$ is far from the observed counts $y_i$, and the chain can take thousands of iterations to escape. Monitor trace plots for each ${\beta }_j$; if the acceptance rate is too low, shrink $c$; if too high, grow it. With a properly-scaled Laplace covariance the chain typically mixes in a few thousand iterations rather than millions, because the proposal shape matches the posterior shape.

Derivation

Goal. Show that the MH transition kernel $P(x,y)$ satisfies detailed balance with respect to $\pi$, i.e., $\pi (x)P(x,y)=\pi (y)P(y,x)$ for all $x,y$. Together with irreducibility of the chain, this establishes $\pi$ as the unique stationary distribution and, via the ergodic theorem, justifies the Monte Carlo approximation $\frac{1}{N}\sum g({\theta }^{(t)})\to \int gd\pi$.

Write the MH transition kernel as
$P(x,y)=Q(x,y)\alpha (x,y)+{\delta }_x(y)r(x),$
where $r(x)=1-\int Q(x,y^{\prime })\alpha (x,y^{\prime })d{y}^{\prime }$ is the total rejection probability at $x$ and ${\delta }_x$ is a point mass. The rejection term trivially satisfies detailed balance because $\pi (x){\delta }_x(y)=\pi (y){\delta }_y(x)$ (both vanish unless $x=y$), so it suffices to check balance for the accept term in the case $y\ne x$:
$\pi (x)Q(x,y)\alpha (x,y)\stackrel{?}{=}\pi (y)Q(y,x)\alpha (y,x).$

Substituting $\alpha (x,y)=\min (1,\pi (y)Q(y,x)/\pi (x)Q(x,y))$ and distributing the $\min$ through the product,
$\pi (x)Q(x,y)\min \left(1,\frac{\pi (y)Q(y,x)}{\pi (x)Q(x,y)}\right)=\min (\pi (x)Q(x,y),\pi (y)Q(y,x)).$
The right-hand side is symmetric under the swap $x\leftrightarrow y$, so it equals $\min (\pi (y)Q(y,x),\pi (x)Q(x,y))$, which by the same algebra equals $\pi (y)Q(y,x)\alpha (y,x)$. This proves detailed balance.

The acceptance probability $\alpha (x,y)=\min (1,\pi (y)Q(y,x)/\pi (x)Q(x,y))$ was not pulled from thin air: it is the maximal function $\alpha :S\times S\to [0,1]$ that makes the balance equation hold. Any $\alpha$ smaller than this minimum-to-one form would also balance the flow, but would reject more often and mix more slowly — the MH choice is the uniquely efficient one.

See Also

• 3.11.2 Gibbs Sampling — sibling MCMC algorithm, the special case of MH with full-conditional proposals always accepted.
• 3.11 MCMC — parent index: chain theory, ergodic theorem, convergence diagnostics.
• 7.4.2. Rejection Sampling — i.i.d. alternative that requires a dominating proposal and a known dominance constant.
• 3.9 Laplace Approximation — standard way to design proposal covariance ${\Sigma }_{\text{prop}}$ for random walk MH.
• 5.4 Reversibility — the detailed-balance condition MH exploits.
• 5.3 Distribution of Markov Chains — stationary distributions and the ergodic theorem.
• 5.1 Markov Chains — foundational prerequisite.