Metropolis-Hastings algorithm

1 Introduction

1.1 Historical context and motivation

The Metropolis–Hastings algorithm was developed to address a recurring problem in statistics and physics: how to sample from a probability distribution when direct drawing is impractical. Early work in computational physics used Markov chains to simulate equilibrium behavior, and later researchers formalized a general framework that applies to arbitrary target distributions provided only up to a normalization constant.

1.2 Where Metropolis–Hastings fits in MCMC

Markov chain Monte Carlo (MCMC) methods generate dependent samples from a target distribution by constructing a Markov chain whose stationary distribution is the target. Metropolis–Hastings is one of the most widely used general-purpose MCMC techniques because it can handle targets specified through unnormalized densities and it allows flexible proposal mechanisms.

1.3 Core idea: proposal, acceptance, and stationarity

At each iteration, the method proposes a new state using a proposal distribution. The proposal is then accepted with a probability designed so that, once the chain has reached its long-run regime, the distribution of states matches the target. The algorithm’s correctness is tied to properties such as detailed balance or reversibility, which guarantee the target as the stationary distribution.

2 Mathematical Foundations

2.1 Markov chains and stationary distributions

A Markov chain is a sequence of random variables where the next state depends only on the current state. If the chain is run for a long time, the distribution of the current state may converge to a stationary distribution, which remains unchanged by further transitions.

2.1.1 Detailed balance and reversibility

A common sufficient condition for stationarity is reversibility. A chain is reversible with respect to a target density π if, for all pairs of states, the probability flow from x to y equals the flow from y to x: \[ \pi(x) q(y\mid x)\alpha(x,y)=\pi(y) q(x\mid y)\alpha(y,x), \] where q is the proposal density and α is the acceptance probability. When this condition holds (and the chain is suitably irreducible), π is stationary.

2.2 Target distribution and unnormalized densities

In many applications, the target distribution π(x) is known only up to a multiplicative constant. Let \[ \pi(x) \propto \tilde{\pi}(x), \] where \(\tilde{\pi}\) is an unnormalized density. Metropolis–Hastings uses ratios of \(\tilde{\pi}\) values, so the unknown constant cancels, enabling sampling without computing normalization.

2.3 Proposal distributions

The proposal distribution q(y	x) defines how candidate moves are generated given the current state x. Proposals can be symmetric (q(y	x)=q(x	y)), such as random-walk proposals, or asymmetric, such as independence proposals that draw candidates from a fixed distribution.

Proposal design affects both correctness (through the required acceptance adjustment) and efficiency (through how often proposals are accepted and how quickly the chain explores).

2.4 Acceptance probability derivation

Metropolis–Hastings selects acceptance probabilities to ensure that the Markov chain has π as its stationary distribution.

2.4.1 General form (asymmetric proposals)

For a current state x and proposed state y, a standard choice of acceptance probability is \[ \alpha(x,y)=\min\left(1,\ \frac{\pi(y)\,q(x\mid y)}{\pi(x)\,q(y\mid x)}\right). \] Using unnormalized densities \(\tilde{\pi}\), this becomes \[ \alpha(x,y)=\min\left(1,\ \frac{\tilde{\pi}(y)\,q(x\mid y)}{\tilde{\pi}(x)\,q(y\mid x)}\right). \] This form compensates for any asymmetry in the proposal mechanism.

2.4.2 Special cases and simplifications

If the proposal is symmetric, q(y	x)=q(x	y), the proposal ratio cancels and the acceptance rule reduces to

\[ \alpha(x,y)=\min\left(1,\ \frac{\pi(y)}{\pi(x)}\right). \] This is the basis of the simpler Metropolis algorithm for symmetric proposals.

3 Algorithm Description

3.1 Step-by-step procedure

Choose an initial state \(x_0\).
For iteration t=0,1,2,…:

Sample a proposal \(y \sim q(\cdot\mid x_t)\).
Compute the acceptance probability

\[ \alpha(x_t,y)=\min\left(1,\ \frac{\tilde{\pi}(y)\,q(x_t\mid y)}{\tilde{\pi}(x_t)\,q(y\mid x_t)}\right). \]

Draw u from Uniform(0,1).
If \(u \le \alpha(x_t,y)\), set \(x_{t+1}=y\); otherwise set \(x_{t+1}=x_t\).

The resulting sequence \(\{x_t\}\) constitutes samples from the chain.

3.2 Choice of initial state

The initial state affects early iterations. If the chain is run long enough and satisfies conditions for convergence, its distribution eventually becomes insensitive to the starting point. In practice, one often discards an initial segment (burn-in) to reduce influence from transient behavior.

3.3 Tuning parameters and practical considerations

Tuning commonly refers to adjusting the scale or shape of the proposal distribution. Too small a proposal can lead to slow exploration (high correlation between successive states), while too large a proposal can reduce acceptance probability, also slowing progress.

A target acceptance rate is sometimes used as a heuristic, but the optimal choice depends on problem dimension, geometry, and the proposal family.

3.4 Burn-in, convergence, and sample dependence

MCMC samples are correlated. A typical workflow includes:

running the chain for enough iterations,
discarding early draws to mitigate burn-in,
assessing whether the remaining samples represent the stationary regime.

Because draws depend on prior states, statistical uncertainty must account for autocorrelation rather than assuming independence.

4 Convergence and Diagnostics

4.1 Mixing behavior and effective sample size

Mixing refers to how rapidly the chain transitions between regions of high probability. Even if the chain is stationary, poor mixing yields highly correlated samples and reduces the number of effectively independent observations.

Effective sample size summarizes this loss by inflating Monte Carlo variance relative to independent sampling.

4.2 Common convergence diagnostics

Diagnostics aim to detect whether the chain has reached its equilibrium distribution and whether results are stable.

4.2.1 Trace plots and autocorrelation

A trace plot shows the sampled state value versus iteration. Stable, stationary-like traces without long drifts are a positive sign, while repeated trends may indicate non-convergence. Autocorrelation functions quantify dependence across lags and help estimate effective sample size.

4.2.2 Multiple chains and comparison

Running multiple chains from dispersed initial points provides evidence of convergence to a common stationary regime. When chains agree in distributional behavior (for example, through visual checks or summary statistics), confidence increases.

4.3 Assessing bias and stability

Even without a definitive proof of convergence in a finite run, practitioners evaluate whether Monte Carlo estimates stabilize as the chain length increases. Sensitivity analyses—such as changing proposal tuning or comparing across chains—can reveal persistent discrepancies that suggest bias.

5 Implementation in Practice

5.1 Selecting proposal families

The proposal family determines how candidates are generated and strongly influences efficiency.

5.1.1 Random-walk proposals

A random-walk proposal proposes \[ y=x+\epsilon, \] where ε is drawn from a distribution centered at zero. This approach is simple but may struggle in complex, high-curvature targets unless step sizes are well matched.

5.1.2 Independence proposals

An independence proposal samples y from a distribution q(y) independent of the current state x. Acceptance then corrects for mismatch between q and the target. This can be efficient when q approximates π well, but it can be ineffective when q places probability mass in unhelpful regions.

5.2 Scaling and acceptance-rate considerations

Practical tuning often adjusts a scale parameter controlling proposal magnitude. Increasing the scale generally decreases acceptance probability but may improve exploration if acceptance remains sufficient. Decreasing the scale increases acceptance but can trap the chain in local moves, producing strong autocorrelation.

5.3 High-dimensional settings

In high dimensions, naive tuning often leads to poor scaling. Proposal covariance matching to the target’s geometry (or adapting proposals) can substantially improve performance. Without such adjustments, the chain may exhibit slow traversal across thin or correlated directions.

5.4 Numerical stability and computation of ratios

Acceptance probabilities depend on ratios involving \(\tilde{\pi}\) and proposal densities. When densities are extremely small, direct computation can underflow. Implementations typically compute in logarithmic space: \[ \log \alpha = \min\left(0,\ \log \tilde{\pi}(y)-\log \tilde{\pi}(x) + \log q(x\mid y)-\log q(y\mid x)\right), \] then exponentiate safely as needed. This improves stability and avoids floating-point errors.

6 Extensions and Variants

The Metropolis algorithm is recovered when the proposal is symmetric. In that setting, the acceptance probability depends only on the ratio π(y)/π(x), simplifying computation while preserving the same stationarity logic.

6.2 Metropolis-within-Gibbs

When the target factors into conditional distributions, Gibbs sampling updates components one at a time. If some conditional updates cannot be sampled directly, Metropolis–Hastings steps can be inserted for those components. The result is a Metropolis-within-Gibbs scheme that combines structured conditional updates with MCMC correction.

6.3 Adaptive Metropolis–Hastings

Adaptive variants modify proposal parameters during sampling based on past states. The adaptation aims to learn a more suitable proposal shape or scale. Correctness typically requires that adaptation diminishes over time, so the chain eventually behaves like a standard Markov chain targeting π.

6.4 Delayed acceptance and approximate proposals

Delayed-acceptance methods use a cheap approximate criterion to reject obviously poor proposals early, saving expensive computations of \(\tilde{\pi}\). Accepted proposals then undergo a second, more accurate acceptance test. This can reduce computational cost when likelihood evaluation is expensive.

6.5 Hamiltonian and other MCMC connections

Metropolis–Hastings provides a template for other MCMC methods. Hamiltonian Monte Carlo (HMC) can be interpreted as using proposals generated from simulated dynamics and then corrected using an acceptance step analogous in spirit to Metropolis–Hastings. While mechanics differ, the common theme is proposal generation followed by acceptance to ensure the target distribution is preserved.

7 Statistical Uses

7.1 Bayesian inference with posterior sampling

In Bayesian analysis, Metropolis–Hastings is commonly used to sample from a posterior distribution proportional to the product of likelihood and prior. Posterior draws enable estimation of parameters and credible intervals when conjugacy is absent or models are too complex for direct sampling.

7.2 Estimating expectations and integrals

For a function h(x), Monte Carlo estimates use sample averages: \[ \mathbb{E}_{\pi}[h(X)] \approx \frac{1}{N}\sum_{i=1}^{N} h(x_i). \] Under appropriate conditions, these estimates converge to the desired expectation as N grows, with uncertainty determined by autocorrelation.

7.3 Posterior predictive and marginalization

Predictive quantities often require integrating over unknown parameters. Samples from the posterior can be propagated through a predictive model to approximate expectations for future observations. Similarly, marginal distributions for subsets of parameters can be obtained by examining the corresponding components of MCMC draws.

7.4 Uncertainty quantification from samples

Credible intervals, tail probabilities, and derived uncertainty measures can be estimated directly from the empirical distribution of the MCMC samples. When dependence is accounted for, the resulting uncertainty reflects both sampling variability and correlation-induced effective sample size reduction.

8 Common Pitfalls

8.1 Poor proposal design and low acceptance rates

If proposals rarely get accepted, the chain remains near its current state for long stretches, leading to slow exploration and biased estimates if run length is insufficient. Conversely, if proposals accept too frequently with tiny moves, the chain may still mix poorly due to strong local correlation.

8.2 Misinterpreting non-convergence

A chain that has not reached stationarity can produce misleading results that look stable over short windows. It is common to mistake early stabilization for convergence. Multiple diagnostics and longer runs are often needed to increase reliability.

8.3 Incorrect use of densities vs. probabilities

Metropolis–Hastings uses density values (or probability mass functions in discrete spaces) in the acceptance ratio. Confusion between densities and probabilities, or ignoring Jacobian factors when changing variables, can lead to incorrect stationary distributions.

8.4 Autocorrelation and misleading effective sample sizes

High autocorrelation can substantially inflate uncertainty. Reporting only the nominal number of samples without accounting for dependence may overstate precision. Effective sample size and uncertainty estimates should reflect the correlation structure.

9 Worked Example (Template)

9.1 Problem setup and target specification

Assume a target distribution for parameter x with density known up to normalization: \[ \pi(x) \propto \tilde{\pi}(x). \] The function \(\tilde{\pi}\) is evaluated for candidate values proposed during sampling. If x is multidimensional, \(\tilde{\pi}\) returns the unnormalized posterior (or other target) density for the full vector.

9.2 Proposal choice and acceptance calculation

Select a proposal distribution q(y	x). For a random-walk normal proposal in one dimension, one might use

\[ y=x+\epsilon,\quad \epsilon\sim\mathcal{N}(0,\sigma^2), \] which is symmetric. The acceptance probability then simplifies to \[ \alpha(x,y)=\min\left(1,\ \frac{\tilde{\pi}(y)}{\tilde{\pi}(x)}\right). \]

For asymmetric proposals, the ratio must include q(x	y)/q(y	x).

9.3 Running the chain and diagnostics

Initialize x0, run for T iterations, discard an initial burn-in period, and then examine trace plots and autocorrelation. If multiple chains are feasible, compare their behavior after burn-in. If diagnostics indicate weak mixing, adjust proposal scale or use a different proposal family and rerun.

9.4 Summarizing results with Monte Carlo estimates

Compute Monte Carlo averages of quantities of interest, such as the posterior mean or functionals h(x). Provide uncertainty intervals derived from the sample distribution, ideally using methods that account for autocorrelation (for example, through effective sample size or batch means).

10.1 Rejection sampling vs. MCMC

Rejection sampling draws independent samples by accepting candidates with probability proportional to the target relative to an envelope distribution. It can be efficient only when the envelope is tight. MCMC trades independence for broader applicability, often requiring careful diagnostics but enabling sampling from difficult targets.

10.2 Importance sampling comparison

Importance sampling reweights samples from a proposal to approximate expectations under the target. When proposals are poorly matched, weights can become highly variable, leading to instability. Metropolis–Hastings similarly relies on a proposal distribution, but instead of reweighting, it corrects via accept/reject steps to create a Markov chain with the correct stationary distribution.

10.3 Other MCMC methods overview

Beyond Metropolis–Hastings, there are many MCMC algorithms, including Gibbs sampling for conditionally tractable models, slice sampling, and Hamiltonian methods for continuous variables. Each method balances ease of implementation, computational cost per iteration, and efficiency in exploration.

10.4 The role of ergodicity in MCMC

Ergodicity links the long-run time averages of a chain to expectations under the stationary distribution. In practice, this property underwrites the use of sample averages as estimators. When ergodicity fails due to reducibility or other structural issues, the chain may explore only part of the state space, preventing correct inference.