Metropolis–Hastings Algorithm

Metropolis–Hastings algorithm is the foundational Markov Chain Monte Carlo (MCMC) method in computation and algorithms for generating a sequence of correlated random samples from a target probability distribution π(x) over a high-dimensional state space when direct sampling or numerical quadrature is intractable. By constructing an ergodic, time-reversible Markov chain whose stationary distribution equals π(x), the algorithm requires evaluating only a function f(x) proportional to the target density—thereby canceling the intractable normalizing partition function Z = ∫ f(x) dx in the acceptance ratio. First developed in 1953 by Nicholas Metropolis, Arianna W. Rosenbluth, Marshall N. Rosenbluth, Augusta H. Teller, and Edward Teller at Los Alamos National Laboratory to simulate hard-sphere fluids under the canonical Boltzmann distribution, and generalized to asymmetric proposal distributions by Canadian statistician W. Keith Hastings in 1970, the algorithm underpins modern Bayesian inference, statistical physics, and statistical phylogenetics.

Eliminating the Intractable Partition Function

In both statistical mechanics and Bayesian inference, the target probability distribution π(x) over a d-dimensional parameter or configuration vector x ∈ Ω ⊆ R^d is known up to a multiplicative normalizing constant Z:

π(x) = f(x) / Z, where Z = ∫_Ω f(x) dx

In canonical statistical mechanics, f(x) = exp(−E(x) / (k_BT)) is the Boltzmann weight of microstate x with potential energy E(x), and Z is the canonical partition function. In Bayesian statistics, f(θ) = P(D | θ)P(θ) is the product of the data likelihood and prior density, while Z = P(D) = ∫ P(D | θ)P(θ) dθ is the marginal likelihood (model evidence). When the dimension d exceeds roughly 10 to 20, deterministic grid quadrature scales as O(M^d) (the curse of dimensionality), and basic rejection or importance Monte Carlo methods fail because independent uniform or prior proposals almost never land inside the narrow high-probability typical set of π(x).

Detailed Balance and the Hastings Acceptance Ratio

The Metropolis–Hastings algorithm circumvents Z by designing a Markov transition kernel P(x → x′) that satisfies the detailed balance (microscopic reversibility) condition with respect to π(x):

π(x) P(x → x′) = π(x′) P(x′ → x), for all x, x′ ∈ Ω

Integrating both sides with respect to x confirms that π(x′) = ∫ π(x)P(x → x′) dx, proving that π is an invariant stationary distribution of the transition kernel. The algorithm separates each transition P(x → x′) (x ≠ x′) into two sub-steps:

1. Proposal Step: Given the current state x_t = x, draw a candidate state x′ from a user-specified conditional proposal distribution q(x′ | x).
2. Acceptance–Rejection Step: Accept the move (x_t+1 = x′) with probability α(x, x′), or reject the move and retain the current state (x_t+1 = x) with probability 1 − α(x, x′).

Substituting P(x → x′) = q(x′ | x) α(x, x′) into the detailed balance equation yields:

α(x, x′) / α(x′, x) = [π(x′) q(x | x′)] / [π(x) q(x′ | x)] = [f(x′) q(x | x′)] / [f(x) q(x′ | x)]

Crucially, the unknown constant Z cancels in the ratio π(x′)/π(x) = f(x′)/f(x). In 1970, W. K. Hastings showed that detailed balance is satisfied while maximizing state mobility by setting the Metropolis–Hastings acceptance probability to:

α(x, x′) = min(1, [f(x′) q(x | x′)] / [f(x) q(x′ | x)])

The factor q(x | x′) / q(x′ | x) is the Hastings ratio, which corrects for any asymmetry in the proposal mechanism. To prevent numerical underflow when likelihoods span hundreds of orders of magnitude, computational implementations evaluate the log-acceptance ratio ln r = [ln f(x′) − ln f(x)] + [ln q(x | x′) − ln q(x′ | x)] and accept whenever ln u < ln r for a uniform draw u ~ Uniform(0, 1). Provided the resulting chain is irreducible (capable of reaching any region of positive probability π > 0) and aperiodic, the ergodic theorem ensures that sample averages (1/N) ∑_t=1^N g(x_t) converge almost surely to the true expectation E_π[g(X)] as N → ∞.

Symmetric Random-Walk Metropolis (1953)

In the original 1953 algorithm of Metropolis, Rosenbluth, Rosenbluth, Teller, and Teller, the proposal distribution is symmetric: q(x′ | x) = q(x | x′), such as an isotropic Gaussian perturbation x′ = x + ε with ε ~ N(0, σ²I). Under symmetry, the Hastings ratio equals 1 and the acceptance probability simplifies to α(x, x′) = min(1, f(x′)/f(x)). In canonical statistical mechanics where f(x) = exp(−E(x)/(k_BT)), any energy-lowering move (ΔE = E(x′) − E(x) ≤ 0) is accepted with probability 1, whereas an uphill energy-raising move (ΔE > 0) is accepted with Boltzmann probability exp(−ΔE/(k_BT)), enabling the sampler to escape local potential energy minima.

For a d-dimensional target distribution, Gareth Roberts, Andrew Gelman, and Walter Gilks (1997) proved via diffusion-limit scaling that as d → ∞, the proposal variance of a Gaussian random-walk Metropolis sampler must scale as σ² ≈ (2.38)²/d, yielding an asymptotic optimal acceptance rate of 0.234 (23.4%; compared to ~0.44 in one dimension).

Gibbs Sampling and Hamiltonian Monte Carlo (HMC)

Two major extensions of the Metropolis–Hastings framework address high-dimensional parameter spaces:

• Gibbs Sampling (Geman & Geman, 1984): Updates one coordinate x_j at a time by drawing from the full conditional distribution q(x′_j | x_j, x_−j) = π(x′_j | x_−j), yielding an acceptance probability α = 1 for every draw.
• Hamiltonian Monte Carlo (HMC; Duane et al., 1987) and the No-U-Turn Sampler (NUTS; Hoffman & Gelman, 2014): Augments x with auxiliary Gaussian momentum variables p ~ N(0, M) and integrates Hamilton's equations using a symplectic leapfrog integrator guided by ∇_x ln f(x). A Metropolis–Hastings acceptance test corrects numerical discretization errors in H(x, p) = −ln f(x) + (1/2)p^TM⁻¹p, converting random-walk diffusion into coherent trajectories across the typical set.

Burn-In, Effective Sample Size (ESS), and Gelman–Rubin R̂

Because MCMC samples are serially correlated and initialized from an arbitrary point x₀, practical inference relies on three diagnostic checks:

1. Burn-In Removal: The initial B warm-up iterations during which the chain relaxes from x₀ into the stationary typical set are discarded before computing posterior expectations.
2. Effective Sample Size (ESS): Under lag-k autocorrelation ρ_k, the sample-mean variance over N MCMC draws equals that of N_eff = N / τ_int independent draws, where τ_int = 1 + 2 ∑_k=1^∞ ρ_k is the integrated autocorrelation time.
3. Gelman–Rubin Potential Scale Reduction Factor (R̂): Using M ≥ 4 chains initialized from overdispersed starting points, Andrew Gelman and Donald Rubin (1992) defined R̂ = √(V̂ / W), comparing pooled variance V̂ to mean within-chain variance W. Values satisfying R̂ < 1.01 (and bulk/tail ESS > 400) confirm convergence across chains.

Bayesian Phylogenetics and Multimodal Tree-Space Exploration

In evolutionary biology, Bayesian phylogenetic suites (MrBayes, BEAST, RevBayes) use the Metropolis–Hastings algorithm to sample the joint posterior over discrete tree topologies τ, continuous branch lengths v, and substitution parameters θ given a sequence alignment D. Because unrooted bifurcating topologies for n taxa grow factorially as (2n − 5)!! (exceeding 2.8 × 10⁷⁶ trees for n = 50), exhaustive enumeration is intractable.

Phylogenetic MCMC chains propose topological rearrangements via Nearest-Neighbor Interchange (NNI), Subtree Pruning and Regrafting (SPR), and Tree Bisection and Reconnection (TBR), evaluating likelihoods P(D | τ, v, θ) via Felsenstein's pruning algorithm and correcting for asymmetric neighborhood sizes through the Hastings ratio. Metropolis-Coupled MCMC (MC³ / parallel tempering) runs a cold chain (β₁ = 1) alongside heated chains sampling flattened power posteriors π(x)^β_k (0 < β_k < 1), proposing periodic Metropolis–Hastings swaps across temperatures to cross multimodal valleys.

Why is an acceptance rate of 100% undesirable in a Random-Walk Metropolis sampler?

An acceptance rate near 100% means the proposal step size σ is too small: every proposed state x′ lies so close to x that f(x′)/f(x) ≈ 1, causing the chain to creep across parameter space via slow Brownian diffusion with high autocorrelation and low Effective Sample Size. Conversely, an acceptance rate near 0% means proposals are too large and step outside the support of π(x). Balancing step size against rejection rate yields the 23.4% optimal acceptance rate in high dimensions.

When is the Hastings correction factor q(x | x′) / q(x′ | x) required?

The Hastings ratio is required whenever the proposal distribution is asymmetric (q(x′ | x) ≠ q(x | x′)). For example, when sampling a positive parameter θ > 0 (such as a branch length or variance) using a log-normal multiplier proposal θ′ = θ exp(ε) with ε ~ Uniform(−δ, δ), the Jacobian of the transformation requires multiplying the acceptance ratio by the Hastings factor θ′/θ to preserve detailed balance.

How does the Metropolis–Hastings algorithm differ from rejection sampling?

Rejection sampling generates independent draws by bounding the target density beneath an envelope M g(x) ≥ f(x), but its acceptance rate decays as O(exp(−c d)) with dimension d. The Metropolis–Hastings algorithm sacrifices sample independence in exchange for local Markovian proposals that track the high-probability manifold across hundreds or thousands of dimensions.

1. Metropolis, N., Rosenbluth, A. W., Rosenbluth, M. N., Teller, A. H., & Teller, E. (1953). Equation of state calculations by fast computing machines. The Journal of Chemical Physics, 21(6), 1087–1092. https://doi.org/10.1063/1.1699114
2. Hastings, W. K. (1970). Monte Carlo sampling methods using Markov chains and their applications. Biometrika, 57(1), 97–109. https://doi.org/10.1093/biomet/57.1.97
3. Gelman, A., & Rubin, D. B. (1992). Inference from iterative simulation using multiple sequences. Statistical Science, 7(4), 457–472. https://doi.org/10.1214/ss/1177011136
4. Roberts, G. O., Gelman, A., & Gilks, W. R. (1997). Weak convergence and optimal scaling of random walk Metropolis algorithms. The Annals of Applied Probability, 7(1), 110–120. https://doi.org/10.1214/aoap/1034625254
5. Ronquist, F., Teslenko, M., van der Mark, P., Ayres, D. L., Darling, A., Höhna, S., Larget, B., Liu, L., Suchard, M. A., & Huelsenbeck, J. P. (2012). MrBayes 3.2: Efficient Bayesian phylogenetic inference and model choice across a large model space. Systematic Biology, 61(3), 539–542. https://doi.org/10.1093/sysbio/sys029