Metropolis-Hastings Algorithm

Learning Objectives

By the end of this section, you will be able to:

Where You'll Apply This: Bayesian inference when conjugacy doesn't exist, sampling from complex posteriors in hierarchical models, uncertainty quantification in ML models, Monte Carlo integration for high-dimensional integrals, and any scenario where direct sampling is impossible.

The Big Picture: Sampling the Unsampleable

In the previous section, we learned that MCMC methods generate samples from probability distributions by constructing a Markov chain whose stationary distribution matches our target. But how do we actually construct such a chain? The Metropolis-Hastings algorithm provides an elegantly simple answer.

The key insight is revolutionary: we don't need to know the normalizing constant of our target distribution. We only need to evaluate the unnormalized density at any point. This is precisely what Bayesian inference requires - we can compute the likelihood times prior, but the evidence (normalizing constant) is often intractable.

Historical Context: From Physics to Statistics

The algorithm has a fascinating origin. In 1953, Nicholas Metropolis and colleagues at Los Alamos National Laboratory developed it to simulate the behavior of liquids on early computers - the MANIAC I. They needed to compute averages over configurations of molecules, but direct computation was impossible due to the astronomical number of states.

In 1970, W.K. Hastings generalized the algorithm to allow asymmetric proposal distributions, vastly expanding its applicability. This generalization is what we call the Metropolis-Hastings algorithm today. The method remained relatively obscure until the 1990s when statisticians recognized its power for Bayesian inference - and it revolutionized computational statistics.

The Algorithm: Step by Step

The Metropolis-Hastings algorithm is remarkably simple. Given a target distribution $\pi(\theta)$ and a proposal distribution $q(\theta' | \theta)$:

The Acceptance Probability

The acceptance probability $\alpha$ has a beautiful interpretation. Let's break it down:

The acceptance probability ensures two crucial properties:

• Moves uphill are always accepted: If $\pi(\theta') > \pi(\theta)$, then $\alpha = 1$
• Moves downhill are sometimes accepted: With probability proportional to the density ratio

Why Does It Work? Detailed Balance

The MH algorithm works because it satisfies detailed balance (also called reversibility):

Here $P(\theta \to \theta') = q(\theta'|\theta) \cdot \alpha(\theta, \theta')$ is the probability of transitioning from θ to θ'. Detailed balance implies that $\pi$ is a stationary distribution of the chain. With additional conditions on the proposal (irreducibility, aperiodicity), the chain converges to $\pi$ from any starting point.

Interactive: Single Step Visualizer

Watch how a single step of Metropolis-Hastings works. The target distribution is a mixture of two Gaussians (bimodal). Click "Step" to see each propose → evaluate → accept/reject cycle.

Interactive: Chain Explorer

Run the sampler continuously and watch the chain build up a histogram that approximates the target distribution. Observe how samples gradually fill in both modes of the bimodal target.

The Proposal Distribution

The choice of proposal distribution $q(\theta'|\theta)$ dramatically affects sampling efficiency. Common choices include:

Tuning the Proposal Width

For random walk proposals, the standard deviation σ creates a fundamental trade-off:

Interactive: Proposal Width Comparison

See the dramatic effect of proposal width on sampling efficiency. Watch three chains run simultaneously with different σ values:

Convergence Diagnostics

How do we know if our chain has converged to the target distribution? Since we can't directly compare samples to the unknown target, we use indirect diagnostics.

Burn-in Period

The chain starts from an arbitrary initial position that may be far from high-probability regions. The burn-in period is the number of initial samples we discard to remove bias from the starting point.

Gelman-Rubin R̂ Statistic

Run multiple chains from different starting points and compare them. The Gelman-Rubin statistic (also called "potential scale reduction factor") measures this:

• R̂ ≈ 1: Chains have converged to the same distribution - good!
• R̂ > 1.1: Chains haven't mixed - run longer or improve proposal
• R̂ >> 1: Serious convergence failure - diagnose and fix

Effective Sample Size

Due to autocorrelation, consecutive MCMC samples are not independent. The effective sample size (ESS) estimates how many truly independent samples your chain is worth:

If your chain has 10,000 samples but ESS = 500, you effectively have only 500 independent samples. High autocorrelation (from small proposal σ) wastes computational effort.

Interactive: Convergence Diagnostics

Watch multiple chains started from very different positions. Initially they're in different regions (high R̂), but as they run, they converge to explore the same distribution (R̂ → 1).

Interactive: Target Distribution Gallery

Different target distributions pose different challenges. Explore how MH performs on various targets, from simple unimodal distributions to challenging multimodal and skewed targets.

Python Implementation

Here's a complete, production-ready implementation of Metropolis-Hastings with examples:

log π(θ|D) ∝ log P(D|θ) + log P(θ)

log(U) < log(α) ⟺ U < α

1import numpy as np
2import matplotlib.pyplot as plt
3
4# =============================================
5# Core Metropolis-Hastings Implementation
6# =============================================
7
8def metropolis_hastings(
9    target_log_prob,    # Log of target density (unnormalized is fine!)
10    proposal_sample,    # Function: current -> proposal
11    x_init,             # Starting point
12    n_samples,          # Number of samples to generate
13    burn_in=1000        # Samples to discard
14):
15    """
16    General Metropolis-Hastings sampler.
17
18    Key insight: We only need log p(x) up to a constant,
19    since we compute ratios p(x')/p(x).
20    """
21    samples = []
22    x_current = x_init
23    acceptances = 0
24
25    for i in range(n_samples + burn_in):
26        # Step 1: Propose a new state
27        x_proposal = proposal_sample(x_current)
28
29        # Step 2: Compute log acceptance ratio
30        # For symmetric proposals: log α = log p(x') - log p(x)
31        log_alpha = target_log_prob(x_proposal) - target_log_prob(x_current)
32
33        # Step 3: Accept/reject
34        log_u = np.log(np.random.random())
35        if log_u < log_alpha:  # Equivalent to u < exp(log_alpha)
36            x_current = x_proposal
37            if i >= burn_in:
38                acceptances += 1
39
40        # Step 4: Record sample (after burn-in)
41        if i >= burn_in:
42            samples.append(x_current)
43
44    acceptance_rate = acceptances / n_samples
45    return np.array(samples), acceptance_rate
46
47
48# =============================================
49# Example 1: Sampling from a Bimodal Distribution
50# =============================================
51
52def bimodal_log_prob(x):
53    """Log probability of mixture of two Gaussians."""
54    mode1, mode2 = -3, 3
55    sigma = 1.0
56    # Log of sum requires log-sum-exp trick for stability
57    log_g1 = -0.5 * ((x - mode1) / sigma) ** 2
58    log_g2 = -0.5 * ((x - mode2) / sigma) ** 2
59    # log(0.5 * exp(log_g1) + 0.5 * exp(log_g2))
60    max_log = max(log_g1, log_g2)
61    return max_log + np.log(0.5 * np.exp(log_g1 - max_log) +
62                            0.5 * np.exp(log_g2 - max_log))
63
64# Symmetric Gaussian proposal
65proposal_sigma = 2.5
66
67def gaussian_proposal(x_current):
68    return x_current + np.random.normal(0, proposal_sigma)
69
70# Run the sampler
71samples, acc_rate = metropolis_hastings(
72    target_log_prob=bimodal_log_prob,
73    proposal_sample=gaussian_proposal,
74    x_init=0.0,
75    n_samples=10000,
76    burn_in=1000
77)
78
79print(f"Acceptance rate: {acc_rate:.1%}")
80print(f"Sample mean: {np.mean(samples):.3f}")
81print(f"Sample std: {np.std(samples):.3f}")
82
83
84# =============================================
85# Example 2: Bayesian Inference for Normal Mean
86# =============================================
87
88def bayesian_normal_mean():
89    """
90    Sample from posterior of Normal mean with known variance.
91
92    Prior: μ ~ N(0, τ²)
93    Likelihood: X_i | μ ~ N(μ, σ²)
94    Posterior: μ | X ~ N(posterior_mean, posterior_var)
95    """
96    # Generate synthetic data
97    np.random.seed(42)
98    true_mu = 2.5
99    sigma = 1.0  # Known likelihood std
100    data = np.random.normal(true_mu, sigma, size=30)
101
102    # Prior parameters
103    prior_mean = 0.0
104    prior_std = 3.0
105
106    # Log posterior (unnormalized)
107    def log_posterior(mu):
108        # Log prior: -0.5 * ((mu - prior_mean) / prior_std)^2
109        log_prior = -0.5 * ((mu - prior_mean) / prior_std) ** 2
110        # Log likelihood: sum of -0.5 * ((x_i - mu) / sigma)^2
111        log_likelihood = -0.5 * np.sum(((data - mu) / sigma) ** 2)
112        return log_prior + log_likelihood
113
114    # Proposal
115    def proposal(mu):
116        return mu + np.random.normal(0, 0.5)
117
118    # Sample
119    samples, acc_rate = metropolis_hastings(
120        target_log_prob=log_posterior,
121        proposal_sample=proposal,
122        x_init=0.0,
123        n_samples=5000,
124        burn_in=500
125    )
126
127    # Compare with analytical posterior
128    n = len(data)
129    post_precision = 1/prior_std**2 + n/sigma**2
130    post_var = 1/post_precision
131    post_mean = post_var * (prior_mean/prior_std**2 + n*np.mean(data)/sigma**2)
132
133    print(f"\nBayesian Normal Mean Inference:")
134    print(f"True μ: {true_mu}")
135    print(f"MH Sample Mean: {np.mean(samples):.4f}")
136    print(f"Analytical Posterior Mean: {post_mean:.4f}")
137    print(f"Acceptance Rate: {acc_rate:.1%}")
138
139bayesian_normal_mean()
140
141
142# =============================================
143# Convergence Diagnostics
144# =============================================
145
146def gelman_rubin(chains):
147    """
148    Compute Gelman-Rubin R-hat statistic for convergence.
149
150    chains: list of arrays, each array is one chain's samples
151    R-hat ≈ 1 indicates convergence; R-hat > 1.1 suggests more sampling needed
152    """
153    m = len(chains)  # Number of chains
154    n = min(len(c) for c in chains)  # Use shortest chain length
155
156    # Chain means
157    chain_means = [np.mean(c[:n]) for c in chains]
158    overall_mean = np.mean(chain_means)
159
160    # Between-chain variance
161    B = n / (m - 1) * sum((cm - overall_mean)**2 for cm in chain_means)
162
163    # Within-chain variance
164    chain_vars = [np.var(c[:n], ddof=1) for c in chains]
165    W = np.mean(chain_vars)
166
167    # Pooled variance estimate
168    var_est = (n - 1) / n * W + B / n
169
170    # R-hat
171    R_hat = np.sqrt(var_est / W) if W > 0 else np.inf
172    return R_hat

Important Variants

AI/ML Applications

Metropolis-Hastings and its variants are foundational to modern probabilistic machine learning:

Common Pitfalls

Knowledge Check

Test your understanding of the Metropolis-Hastings algorithm with this interactive quiz.

Summary

Looking Ahead: In the next section, we'll explore Gibbs Sampling - a special case of MH that uses conditional distributions as proposals. When full conditionals are available, Gibbs sampling provides a systematic way to sample from complex joint distributions one coordinate at a time.