Metropolis-Hastings Algorithm

Definition

The Metropolis-Hastings (MH) algorithm is a general-purpose Markov Chain Monte Carlo (MCMC) method that generates samples from a target distribution $\pi(\theta)$ by iterating a proposal–accept/reject step: at each iteration draw a candidate $\theta^*$ from a proposal density $q(\theta^*|\theta^{t-1})$ , then accept it with probability

$r = \min\!\left(1,\; \frac{\pi(\theta^*)\, q(\theta^{t-1} \mid \theta^*)}{\pi(\theta^{t-1})\, q(\theta^* \mid \theta^{t-1})}\right)$

If rejected, the chain stays at $\theta^{t-1}$ . The ratio corrects for asymmetric proposals so that the chain's invariant distribution is exactly $\pi$ , regardless of the proposal form. The acceptance function is derived from the reversibility condition (detailed balance): requiring $\pi(x)p(x,y) = \pi(y)p(y,x)$ and minimizing rejection probability yields the Hastings form uniquely (Chib-Greenberg 1995; Peskun 1973).

Key Ideas

The acceptance ratio $r$ requires only the ratio $\pi(\theta^*)/\pi(\theta^{t-1})$ , so the normalising constant of $\pi$ cancels — this is the key advantage for Bayesian posteriors.
The Gibbs Sampler is a special case: when $q(\theta^*|\theta^{t-1}) = \pi(\theta^* \mid \theta^{t-1}_{-i}, Y)$ (the full conditional), the ratio $r = 1$ everywhere and every candidate is accepted.
Peskun (1973) optimality. Among all acceptance functions $\alpha(x,y)$ that leave $\pi$ invariant, the Hastings form is optimal: it minimises the rejection probability for any pair $(x,y)$ and therefore minimises the asymptotic variance $\sigma^2$ of MCMC mean estimates. Barker's (1965) form $\pi(y)q(y,x)/[\pi(x)q(x,y) + \pi(y)q(y,x)]$ is strictly worse.
Product of Kernels principle (Chib-Greenberg 1995). If kernel $P_1$ preserves $\pi(\theta_1 \mid \theta_2)$ and $P_2$ preserves $\pi(\theta_2 \mid \theta_1)$ , then the composition $P_1 P_2$ preserves the joint $\pi(\theta_1, \theta_2)$ . This justifies all Metropolis-within-Gibbs hybrids: blocks with conjugate conditionals use $\alpha = 1$ (pure Gibbs); intractable blocks use an MH step — the composition still targets $\pi$ .
Convergence and a Central Limit Theorem (CLT) follow from general state-space Markov chain theory: irreducibility and aperiodicity guarantee ergodic convergence; geometric ergodicity (verified via a drift condition) validates confidence intervals.

How It Works

The Algorithm

At iteration $t$ :

Draw $\theta^* \sim q(\cdot \mid \theta^{t-1})$ .
Compute $r = \min(1,\; \pi(\theta^*)q(\theta^{t-1}|\theta^*) / [\pi(\theta^{t-1})q(\theta^*|\theta^{t-1})])$ .
Set $\theta^t = \theta^*$ with probability $r$ ; else $\theta^t = \theta^{t-1}$ .

Discard an initial burn-in; use $\{\theta^t\}$ as approximate draws from $\pi$ .

Derivation from Reversibility (Chib-Greenberg 1995)

The MH acceptance function follows from the detailed balance (reversibility) condition: the chain's transition kernel $p(x,y)$ must satisfy

$\pi(x)\,p(x,y) = \pi(y)\,p(y,x) \quad \forall\, x, y$

For a proposal $q(\cdot|\cdot)$ , write $p(x,y) = q(x,y)\,\alpha(x,y)$ for $y \neq x$ . Substituting and rearranging:

$\pi(x)\,q(x,y)\,\alpha(x,y) = \pi(y)\,q(y,x)\,\alpha(y,x)$

The Peskun-optimal $\alpha$ that satisfies this and maximises $\alpha(x,y) + \alpha(y,x)$ (minimising total rejection) is

$\alpha(x,y) = \min\!\left(1,\, \frac{\pi(y)\,q(y,x)}{\pi(x)\,q(x,y)}\right)$

When $\pi(y)q(y,x) > \pi(x)q(x,y)$ , set $\alpha(x,y) = 1$ and $\alpha(y,x) = \pi(x)q(x,y)/[\pi(y)q(y,x)] < 1$ .

Five Candidate-Generating Families (Chib-Greenberg 1995)

Family	Proposal	Acceptance ratio
Random walk	$y = x + \varepsilon$ , $q$ symmetric	$\min\{1,\, \pi(y)/\pi(x)\}$
Independence chain	$q(y \mid x) = q(y)$ , independent of $x$	$\min\{1,\, \pi(y)q(x)/[\pi(x)q(y)]\}$
Rejection-sampling	$\pi \propto \psi h$ , propose from $h$	$\min\{1,\, \psi(y)/\psi(x)\}$ ; Gibbs = limit
Autoregressive / reflection	$y = a + B(x-a) + z$ ; $B = -I$ reflects	Full Hastings ratio; negative serial corr.
Block-at-a-time	Partition $\theta = (\theta_1, \theta_2, \ldots)$ ; update blocks in turn	Block-specific Hastings ratios

M-H Accept-Reject (A-R) Algorithm (Chib-Greenberg 1995, §6.1)

When the target $f(\cdot)$ is bounded above by $c \cdot h(\cdot)$ on a set $C = \{x: f(x) < c\,h(x)\}$ , a pseudodominating density

$q(y) = \begin{cases} f(y)/(cd) & y \in C \\ h(y)/d & y \notin C \end{cases}$

( $d$ chosen so $q$ integrates to 1) can be used as the MH proposal. The acceptance probabilities for four cases $(x \in/\notin C, y \in/\notin C)$ are tabulated in Chib-Greenberg (1995) Table 1. The algorithm combines the efficiency of rejection sampling inside the "easy" region $C$ with a Metropolis correction outside it, ensuring invariance of $\pi$ without burn-in bias.

Six Metropolis Kernel Types (Tierney 1994)

Tierney (1994) unifies the MH literature under six kernel types:

Kernel	Proposal $q(\cdot \mid \theta)$	Notes
Random-walk	$\theta^* = \theta + \varepsilon$ , $\varepsilon \sim \mathcal{N}(0,\tau^2)$	Most common; tune $\tau$ for acceptance rate $0.44$ (scalar) / $0.23$ (high-dim)
Independence chain	$q(\theta^* \mid \theta) = q(\theta^*)$ — independent of current state	Efficient when $q \approx \pi$ ; collapses when tails of $q$ are too thin
Rejection-sampling chain	Embeds a reject-sampling step in MH; guarantees exactly $\pi$ even when the envelope occasionally fails	Useful for bounded supports
Autoregressive / shrinkage	$\theta^* = \rho\,\theta + \varepsilon$ , $\varepsilon \sim \mathcal{N}(0,\sigma^2_\varepsilon)$	Shrinks proposals toward origin; useful for near-stationary targets
Grid-based (Griddy)	Evaluate $\pi$ on a grid; cumulative distribution function (CDF) inversion to propose	Embeds Griddy-Gibbs within MH framework
Reflection	$\theta^* = 2a - \theta$ reflected around $a$	Induces negative serial correlations, reducing asymptotic variance

Convergence Theory (Tierney 1994)

Drawing on Nummelin (1984) general state-space theory:

Theorem 1 (Ergodic convergence). If the chain is $\pi$ -irreducible and aperiodic, then $\pi$ is the unique invariant distribution and for $\pi$ -a.e. starting point $x$ : $\|P^n(x,\cdot) - \pi(\cdot)\| \to 0$

Geometric ergodicity via drift condition. If there exist $V: E \to [1,\infty)$ , a "small" set $C$ , and constants $\lambda < 1$ , $b < \infty$ such that: $PV(x) \leq \lambda V(x) + b\cdot\mathbf{1}_C(x)$ then $\|P^n(x,\cdot) - \pi\| \leq M(x)\rho^n$ for some $\rho < 1$ . This is the standard sufficient condition for CLT validity.

Theorem 4 (CLT). If the chain is geometrically ergodic and $g \in L^{2+\varepsilon}(\pi)$ : $\sqrt{n}\!\left(n^{-1}\sum_{j=1}^n g(\theta^j) - E_\pi[g]\right) \xrightarrow{d} \mathcal{N}(0,\sigma^2)$

The asymptotic variance is the spectral density of the chain at zero frequency: $\sigma^2 = \mathrm{Var}_\pi(g) + 2\sum_{j=1}^\infty \mathrm{Cov}_\pi(g(\theta^0),\, g(\theta^j))$

Geometric ergodicity guarantees this sum converges absolutely. Uniform ergodicity (stronger) requires only $g \in L^2(\pi)$ (Theorem 5).

Metropolis-within-Gibbs

When a Gibbs sampler has one or more blocks with non-standard full conditionals, MH steps can be nested inside the Gibbs cycle. Each block that cannot be drawn directly is replaced by one or more MH steps targeting that block's full conditional; blocks with standard conjugate forms continue as direct draws. The resulting algorithm is a hybrid (mixture/cycle of kernels) that retains $\pi$ as its invariant distribution. See Time-Varying Parameter VAR for the five-block Metropolis-within-Gibbs sampler of Cogley-Morozov-Sargent (2003) and Noninformative Prior for VAR for the hit-and-run variant.

WLS Proposal for GLMMs (Gamerman 1997)

For generalized linear mixed models (GLMMs), a highly efficient weighted least squares (WLS) Metropolis proposal is obtained from the iteratively reweighted least squares (IWLS) linearization of the non-Gaussian likelihood. At state $\eta^{(t)}$ , compute pseudo-response $\tilde{y}_i = \eta_i + (y_i - \mu_i)g'(\mu_i)$ and weight $W_i^{-1} = b''(\theta_i)[g'(\mu_i)]^2$ . The resulting working-model posterior $N(m^{(t)}, C^{(t)})$ — a generalized least squares (GLS) posterior for the linearized Normal problem — serves as an independence-chain proposal whose acceptance rate exceeds 85–98% in practice. This is an instance of the independence chain family where the proposal is locally calibrated to the posterior curvature rather than fixed. See Generalized Linear Mixed Model for full details.

Prefetching: Parallelising MH via a Metropolis Tree (Strid 2010)

Standard MH is intrinsically serial: each evaluation of $\pi(\theta^*)$ must complete before the next proposal can be formed. Prefetching breaks this bottleneck by speculating on the accept/reject outcome and evaluating the posterior at multiple future candidate states in parallel.

Metropolis tree. With $P$ processors, build a binary tree where each node corresponds to a possible future state: the left branch (acceptance) moves to $\theta^*$ ; the right branch (rejection) stays at $\theta^{t-1}$ . A tour of size $P$ selects $P$ nodes and evaluates them in parallel. The expected draws per tour: $D(P) = \sum_{p=1}^P \Pr(\text{node } p \text{ selected})$ is the theoretical speedup over one serial step.

Five variants. Algorithms differ in how branch probabilities are assigned (which determines $D$ ):

Variant	Branch probabilities	$D$ at $P=7$ , $\alpha=0.25$
Basic (Brockwell 2006)	All 0.5	3.00; lower bound $\log_2(P+1)$
Static	Known $\alpha$ : reject = $1-\alpha$ , accept = $\alpha$	3.54
Dynamic (uniform)	Empirical $\alpha_k$ from pilot histogram over $U[0,1]$ bins	3.81
Dynamic (posterior approx.)	$\min(\beta,\, p^(\theta_\text{acc})/p^(\theta_\text{cur}))$ with $p^* = \mathcal{N}(\theta_m, H_m^{-1})$	4.58
Most likely path	Run $P$ steps ahead with $p^*$ ; non-branching tour	4.99

Optimal acceptance rate shifts downward. The joint expected output $E(\alpha,P) = E^1(\alpha) \times D(\alpha,P)$ is maximised at lower $\alpha$ than the serial optimum 0.234 because lower $\alpha$ lengthens accepted-state runs, increasing $D$ :

$P$	$\alpha_\text{opt,P}$
1	0.234
7	0.120
$\infty$	0

Smets-Wouters benchmark (Table 1, $\alpha=0.25$ ). Most likely path achieves theoretical speedup $4.99\times$ ( $P=7$ ) and actual cluster speedup $5.15\times$ on a Lenngren high-performance computing (HPC) cluster with 8 processors (80 min serial → 15 min parallel). Independence chain MH is embarrassingly parallel ( $S(P)\approx P$ ) but achieves only 13% acceptance on the 40-parameter Smets-Wouters (SW) model, making its parallelism advantage illusory.

Limitations. Practical ceiling $\approx 10$ – $15$ processors; most likely path inapplicable with discontinuous likelihoods (nonlinear Dynamic Stochastic General Equilibrium (DSGE)/particle filter); Gibbs steps ( $\alpha=1$ ) incompatible with prefetching.

Tuning and Diagnostics

Acceptance rate target (Roberts-Gelman-Gilks 1994): $\approx 0.45$ for one-dimensional $\theta$ ; $\approx 0.25$ for $d \geq 6$ ; $\approx 0.23$ asymptotically as $d \to \infty$ . These apply to random-walk Gaussian proposals; rates far from target signal a mis-scaled proposal variance.
Burn-in: Discard initial draws until the chain appears stationary (visual inspection of traces; Gelman-Rubin $\hat{R}$ across multiple chains).
Effective sample size: Low ESS relative to total draws indicates slow mixing — reduce proposal variance (random-walk) or improve proposal approximation (independence chain).
Variance estimation: Use the spectral variance formula above; estimate $\sigma^2$ via batch means or a smoothed periodogram (Geweke 1992). See Gibbs Sampler for implementation details.

Why It Matters

MH is the most general MCMC workhorse: it applies whenever the target $\pi$ can be evaluated up to a normalising constant, with no conjugacy requirements.
The Gibbs sampler's theoretical properties (convergence, CLT) are understood as consequences of MH theory — Gibbs is the limit of MH as the proposal becomes the full conditional.
Peskun optimality gives a theoretical justification for the Hastings acceptance form over alternatives, settling a design question that had been open since Barker (1965).
Hybrid MH-Gibbs algorithms power modern time-varying parameter VAR (TVP-VAR), stochastic volatility (SV), and other high-dimensional state-space estimators that mix conjugate and non-conjugate blocks.

Open Questions

Optimal proposal scaling in high dimensions (beyond the $0.23$ rule for isotropic Gaussian proposals) is an active research area; adaptive MCMC methods (Haario-Saksman-Tamminen) auto-tune but require careful justification of ergodicity.
For multimodal targets, standard random-walk MH may mix across modes only rarely; parallel tempering and simulated tempering extend MH to address this but add tuning complexity.
The drift-condition approach verifies geometric ergodicity theoretically; in practice for specific models (e.g. heavy-tailed posteriors), checking the condition can be non-trivial.