Markov Chain Monte Carlo

Introduction

Markov Chain Monte Carlo (MCMC) refers to a family of computational algorithms designed to approximate complex probability distributions by generating a sequence of random samples. Rather than computing intractable integrals analytically, MCMC methods construct a Markov chain whose stationary distribution matches the target distribution of interest. As the chain evolves, the sequence of states converges to the desired distribution, enabling estimation of expectations, variances, and other statistical quantities.

The technique has become foundational in Bayesian statistics, computational physics, machine learning, and computational biology, where it handles high-dimensional parameter spaces that defy direct sampling or grid-based numerical integration.

Historical Development

The origins of MCMC trace back to the Metropolis algorithm (1953), developed to study equilibrium properties of atomic systems in statistical mechanics. Nicholas Metropolis, Arianna Rosenthal, Marshall Rosenbluth, Augusta Teller, and Edward Teller introduced a procedure that accepted or proposed state transitions based on an acceptance probability, ensuring detailed balance with respect to a target Boltzmann distribution.

In 1970, W. K. Hastings generalized the method to asymmetric proposal distributions, yielding the widely used Metropolis–Hastings algorithm. The 1980s saw the introduction of Gibbs sampling by Geman & Geman (1984) and later formalized by Gelfand & Smith (1990), which exploits conditional distributions to simplify sampling in multivariate settings. The 1990s–2000s witnessed rapid expansion into Bayesian computation, with innovations such as Hamiltonian Monte Carlo (HMC) and the No-U-Turn Sampler (NUTS) revolutionizing high-dimensional inference.

Mathematical Framework

Let \(\pi(\mathbf{x})\) denote the target probability density or mass function over a state space \(\mathcal{X}\). In Bayesian contexts, \(\pi(\mathbf{x})\) is typically a posterior distribution proportional to the product of a likelihood and a prior. Direct sampling from \(\pi(\mathbf{x})\) is often impossible due to an intractable normalization constant \(Z\).

MCMC bypasses the need to compute \(Z\) by constructing a Markov chain \(\{\mathbf{x}^{(t)}\}_{t=0}^T\) with transition kernel \(K(\mathbf{x} \to \mathbf{y})\) that satisfies detailed balance (reversibility):

When detailed balance holds, \(\pi(\mathbf{x})\) is the stationary distribution of the chain. Under mild ergodicity conditions, the empirical distribution of the chain converges to \(\pi(\mathbf{x})\) as \(T \to \infty\), allowing Monte Carlo estimation:

Core Algorithms

Metropolis–Hastings (MH)

The MH algorithm proposes a candidate \(\mathbf{y}\) from a proposal distribution \(q(\mathbf{y} | \mathbf{x}^{(t)})\) and accepts it with probability:

If rejected, the chain remains at \(\mathbf{x}^{(t+1)} = \mathbf{x}^{(t)}\). MH is highly flexible but can suffer from slow mixing in high dimensions due to random-walk behavior.

Gibbs Sampling

Gibbs sampling is a special case of MH where the proposal distribution equals the full conditional distribution. For a parameter vector \(\mathbf{x} = (x_1, \dots, x_d)\), each component is updated sequentially:

Gibbs sampling guarantees acceptance probability 1 and excels when conditionals are known and easy to sample from, though it can exhibit strong autocorrelation in highly dependent parameter spaces.

Hamiltonian Monte Carlo (HMC)

HMC leverages concepts from classical mechanics to propose distant states with high acceptance rates. It augments the parameter space with auxiliary momentum variables \(\mathbf{p}\) and simulates Hamiltonian dynamics using the potential energy \(U(\mathbf{x}) = -\log \pi(\mathbf{x})\). The trajectory preserves volume and energy, enabling efficient exploration of correlated distributions. The No-U-Turn Sampler (NUTS) automatically tunes trajectory lengths, making it the default algorithm in modern probabilistic programming systems.

Applications

• Bayesian Inference: Posterior sampling for hierarchical models, latent variable models, and non-conjugate priors.
• Machine Learning: Training Boltzmann machines, variational inference approximation, and Bayesian neural networks.
• Computational Physics: Lattice QCD simulations, protein folding, and statistical mechanics of many-body systems.
• Computational Biology: Phylogenetic tree reconstruction, molecular dynamics, and gene regulatory network inference.
• Operations Research: Stochastic optimization, reliability analysis, and risk assessment in complex systems.

Advantages & Limitations

Strengths: Avoids curse of dimensionality in integral approximation, requires only unnormalized densities, highly adaptable to complex models.

Weaknesses: No guaranteed finite-time convergence, computationally intensive for large datasets, sensitive to proposal tuning, challenging convergence validation in high dimensions.