Markov Chain Monte Carlo

Learning Objectives

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

Where You'll Apply This: Bayesian inference for complex posteriors, generative models, uncertainty quantification in neural networks, probabilistic programming (PyMC, Stan, NumPyro), molecular simulations, and any scenario where direct sampling is intractable.

The Big Picture: Sampling the Impossible

In the previous section on Monte Carlo integration, we learned that we can estimate expectations by averaging over random samples. But there's a catch: we assumed we could directly generate samples from our target distribution. What happens when we can't?

Consider a Bayesian posterior: $P(\theta | D) = \frac{P(D|\theta) P(\theta)}{P(D)}$. The normalizing constant $P(D) = \int P(D|\theta) P(\theta) d\theta$ is typically intractable for complex models. We can evaluate the posterior up to this constant, but we can't sample from it directly.

This is the brilliant insight of Markov Chain Monte Carlo (MCMC): convert a sampling problem into a simulation problem. Instead of sampling directly, we simulate a random process that eventually produces samples from the desired distribution.

Historical Context

MCMC was invented in 1953 by Metropolis, Rosenbluth, Rosenbluth, Teller, and Teller at Los Alamos National Laboratory, originally to simulate molecular systems. Hastings generalized the algorithm in 1970, and the method revolutionized statistics when applied to Bayesian inference in the 1990s with the work of Gelfand and Smith on Gibbs sampling.

Markov Chains: The Foundation

Before diving into MCMC, we need to understand Markov chains - the engine that powers the method.

The Markov Property

A sequence of random variables $X_0, X_1, X_2, \ldots$ is a Markov chainif the future depends only on the present, not the past:

Transition Kernels

The dynamics of a Markov chain are specified by a transition kernel$K(x, x')$ that gives the probability (or density) of moving from state$x$ to state $x'$:

Given a distribution $\mu_n(x)$ at step $n$, the distribution at step $n+1$ is:

Stationary Distributions

A distribution $\pi$ is stationary for a Markov chain if applying one step of the chain leaves the distribution unchanged:

Under mild conditions (irreducibility and aperiodicity), a Markov chain converges to its stationary distribution regardless of where it starts. This is the ergodic theorem - the foundation of MCMC.

Interactive: Convergence to Stationarity

Watch how a probability distribution evolves under repeated application of a Markov chain's transition matrix. Regardless of the starting distribution, the chain converges to its unique stationary distribution.

The Core MCMC Idea

The goal of MCMC is to construct a Markov chain whose stationary distribution equals our target distribution $\pi$ (e.g., the posterior). If we can do this, then running the chain generates samples from $\pi$.

Detailed Balance

A sufficient condition for $\pi$ to be stationary is detailed balance:

Detailed balance implies stationarity because integrating over $x$:

Ergodicity and Convergence Guarantees

For MCMC to work, the chain must be ergodic, meaning it will eventually explore the entire state space. This requires:

With these properties, the ergodic theorem guarantees:

Interactive: MCMC Sampling Visualizer

Watch MCMC in action! This visualization shows a Metropolis-Hastings sampler exploring a bimodal target distribution. Adjust the proposal standard deviation to see how it affects mixing and acceptance rates.

Practical MCMC: Making It Work

Running an MCMC algorithm is only half the battle. Getting reliable samples requires careful attention to several practical issues.

Burn-in Period

Early samples are influenced by the starting point and may not represent the stationary distribution. The burn-in (or warm-up) period allows the chain to "forget" its initialization.

Thinning

Consecutive MCMC samples are correlated. Thinning keeps every $k^{th}$sample to reduce autocorrelation and storage requirements.

Multiple Chains

Running multiple chains from dispersed starting points provides several benefits:

• Convergence diagnostics: Compare chains using R-hat to detect non-convergence
• Mode discovery: Different chains may find different modes in multimodal distributions
• Parallelization: Chains can run independently on multiple cores
• Robustness: If one chain gets stuck, others may still explore well

Interactive: Trace Plots and Diagnostics

Trace plots are the primary visual diagnostic for MCMC convergence. Watch multiple chains explore the parameter space and observe how the histogram converges to the target distribution.

Convergence Diagnostics

"Has my chain converged?" is the fundamental question of MCMC practice. We can never prove convergence, but several diagnostics can detect non-convergence.

R-hat (Gelman-Rubin Diagnostic)

The R-hat statistic compares variance between chains to variance within chains. If chains have converged to the same distribution, these should be similar.

Effective Sample Size

Due to autocorrelation, $N$ MCMC samples provide less information than$N$ independent samples. The Effective Sample Size (ESS)quantifies how many independent samples your correlated chain is worth.

Interactive: Diagnostics Explorer

Explore how different MCMC settings affect convergence diagnostics. Adjust the convergence quality, number of chains, and chain length to see their effects on R-hat and ESS.

Python Implementation

Here's a complete implementation of a basic MCMC sampler with convergence diagnostics. This code demonstrates the core concepts and provides a foundation for more sophisticated samplers.

log π(x) ∝ log P(data|x) + log P(x)

q(x' | x) = N(x'; x, σ²I)

1import numpy as np
2from typing import Callable, Tuple, List
3import matplotlib.pyplot as plt
4
5# ============================================
6# Simple MCMC Framework
7# ============================================
8
9def mcmc_sampler(
10    log_prob: Callable[[np.ndarray], float],
11    initial_state: np.ndarray,
12    n_samples: int,
13    proposal_std: float = 1.0,
14) -> Tuple[np.ndarray, float]:
15    """
16    Generic MCMC sampler using random walk proposals.
17
18    Parameters:
19    -----------
20    log_prob : function
21        Log probability of target distribution (up to constant)
22    initial_state : np.ndarray
23        Starting point for the chain
24    n_samples : int
25        Number of samples to generate
26    proposal_std : float
27        Standard deviation of Gaussian proposal
28
29    Returns:
30    --------
31    samples : np.ndarray
32        Array of shape (n_samples, dim) containing samples
33    acceptance_rate : float
34        Fraction of proposals that were accepted
35    """
36    dim = len(initial_state)
37    samples = np.zeros((n_samples, dim))
38    current = initial_state.copy()
39    current_log_prob = log_prob(current)
40
41    n_accepted = 0
42
43    for i in range(n_samples):
44        # Generate proposal from symmetric Gaussian
45        proposal = current + np.random.normal(0, proposal_std, dim)
46        proposal_log_prob = log_prob(proposal)
47
48        # Compute acceptance probability (Metropolis ratio)
49        log_alpha = proposal_log_prob - current_log_prob
50
51        # Accept or reject
52        if np.log(np.random.uniform()) < log_alpha:
53            current = proposal
54            current_log_prob = proposal_log_prob
55            n_accepted += 1
56
57        samples[i] = current
58
59    acceptance_rate = n_accepted / n_samples
60    return samples, acceptance_rate
61
62# ============================================
63# Example: Sampling from a 2D Mixture of Gaussians
64# ============================================
65
66def bimodal_log_prob(x: np.ndarray) -> float:
67    """Log probability of mixture of two 2D Gaussians."""
68    mu1, mu2 = np.array([-2, 0]), np.array([2, 0])
69    sigma = 1.0
70
71    log_p1 = -0.5 * np.sum((x - mu1)**2) / sigma**2
72    log_p2 = -0.5 * np.sum((x - mu2)**2) / sigma**2
73
74    # Log-sum-exp trick for numerical stability
75    max_log_p = max(log_p1, log_p2)
76    return max_log_p + np.log(
77        np.exp(log_p1 - max_log_p) + np.exp(log_p2 - max_log_p)
78    ) - np.log(2)
79
80# Run MCMC with different proposal scales
81initial = np.array([0.0, 0.0])
82n_samples = 5000
83
84# Small proposal - high acceptance, slow mixing
85samples_small, acc_small = mcmc_sampler(
86    bimodal_log_prob, initial, n_samples, proposal_std=0.3
87)
88print(f"Small proposal: acceptance = {acc_small:.2%}")
89
90# Optimal proposal - balanced mixing
91samples_opt, acc_opt = mcmc_sampler(
92    bimodal_log_prob, initial, n_samples, proposal_std=1.5
93)
94print(f"Optimal proposal: acceptance = {acc_opt:.2%}")
95
96# Large proposal - low acceptance, poor mixing
97samples_large, acc_large = mcmc_sampler(
98    bimodal_log_prob, initial, n_samples, proposal_std=5.0
99)
100print(f"Large proposal: acceptance = {acc_large:.2%}")
101
102# ============================================
103# Convergence Diagnostics
104# ============================================
105
106def effective_sample_size(samples: np.ndarray) -> float:
107    """Compute effective sample size accounting for autocorrelation."""
108    n = len(samples)
109    if n < 10:
110        return n
111
112    mean = np.mean(samples)
113    variance = np.var(samples)
114    if variance == 0:
115        return n
116
117    # Compute autocorrelation
118    max_lag = min(100, n // 2)
119    acf = np.correlate(samples - mean, samples - mean, mode='full')
120    acf = acf[n-1:n-1+max_lag] / (variance * n)
121
122    # Integrated autocorrelation time
123    tau = 1 + 2 * np.sum(acf[1:])
124
125    return n / tau
126
127def gelman_rubin(chains: List[np.ndarray]) -> float:
128    """
129    Compute R-hat (Gelman-Rubin diagnostic).
130
131    R-hat ≈ 1 indicates convergence.
132    R-hat > 1.1 suggests chains haven't converged.
133    """
134    m = len(chains)  # Number of chains
135    n = min(len(c) for c in chains)  # Chain length
136
137    # Chain means and overall mean
138    chain_means = [np.mean(c[:n]) for c in chains]
139    overall_mean = np.mean(chain_means)
140
141    # Between-chain variance
142    B = n / (m - 1) * sum((cm - overall_mean)**2 for cm in chain_means)
143
144    # Within-chain variance
145    chain_vars = [np.var(c[:n], ddof=1) for c in chains]
146    W = np.mean(chain_vars)
147
148    # Pooled variance estimate
149    var_hat = (n - 1) / n * W + B / n
150
151    # R-hat
152    return np.sqrt(var_hat / W)
153
154# Run multiple chains for R-hat
155chains = []
156for start_x in [-5, 0, 5]:
157    initial = np.array([start_x, 0.0])
158    samples, _ = mcmc_sampler(bimodal_log_prob, initial, 2000)
159    chains.append(samples[:, 0])  # First dimension
160
161rhat = gelman_rubin(chains)
162print(f"R-hat: {rhat:.4f}")
163
164# Post burn-in ESS
165burn_in = 500
166ess = effective_sample_size(samples_opt[burn_in:, 0])
167print(f"Effective Sample Size: {ess:.0f} / {n_samples - burn_in}")

Real-World Examples

AI/ML Applications

MCMC and its descendants are fundamental to modern machine learning, especially for uncertainty quantification and generative modeling.

Common Pitfalls

Knowledge Check

Test your understanding of Markov Chain Monte Carlo with this interactive quiz.

Summary

Looking Ahead: In the next section, we'll explore the Metropolis-Hastings algorithm, which provides a general recipe for constructing MCMC transition kernels. This is the foundation of most practical MCMC methods, including the advanced algorithms used in modern probabilistic programming.