Part II · Monte Carlo and MCMC — computing with randomness

draft

If a Markov chain is a machine that walks, Monte Carlo is the discovery that a well-chosen random walk can compute. This chapter follows the idea from a convalescing mathematician's card game to the engine that samples nearly every modern probabilistic model. The theory of the walk itself is chapter 01, Markov chains; here the walk is put to work.

1. The idea: sample instead of solve

Some quantities are defined by a sum or an integral too large to evaluate directly — a high-dimensional integral, the probability of a rare configuration, the expected value of a function over an enormous state space. The Monte Carlo insight is to stop trying to solve the integral and instead estimate it by random sampling: draw many samples, average the thing you care about, and let the law of large numbers (see ch. 01, §1) carry the average to the true value.

To estimate an integral \(I = \int f(x)\,dx\) over a domain, draw \(N\) independent samples and average \(f\) (weighted by the domain measure or sampling density). The sample mean is an unbiased estimator of \(I\). The canonical toy version: to estimate \(\pi\), throw points uniformly at a unit square; the fraction landing inside the inscribed quarter-circle (\(x^2+y^2\le 1\)) approaches \(\pi/4\).

Why it wins: the error is dimension-blind. The standard error of a Monte Carlo estimate scales as \(\sigma/\sqrt{N}\) — error \(\propto 1/\sqrt{N}\), so quadrupling the samples halves the error. Crucially, this rate does not depend on the dimension of the space being integrated. Grid-based quadrature suffers the curse of dimensionality (error \(\sim N^{-k/d}\) in \(d\) dimensions, hopeless for large \(d\)); Monte Carlo pays the same \(1/\sqrt{N}\) whether the space has three dimensions or three million. That single fact is why it dominates statistical physics, finance, and Bayesian inference — all high-dimensional integration in disguise.

2. Origin: solitaire, neutrons, and the bomb

The method was born inside the Manhattan Project. In 1946, recovering from an illness (encephalitis) and playing Canfield solitaire, Stanisław Ulam asked a simple question: what is the probability that a laid-out game comes out? The combinatorics were hopeless — but he realized he could just *deal the game many times and count the wins*. Estimate the probability by repeated random trials.

Ulam described the idea to John von Neumann, and they aimed it at the real problem in front of them: neutron diffusion and transport in fissile material — a high-dimensional integro-differential problem that resisted every deterministic attack. In March 1947 von Neumann wrote to Robert Richtmyer at Los Alamos with a detailed scheme for tracking simulated neutron histories (position, velocity, collision, absorption, fission) by random sampling, explicitly designed for the ENIAC. The first automated fission-core runs followed in 1947–1948.

The name is a code name: Nicholas Metropolis proposed "Monte Carlo" after the Monaco casino — an allusion to an uncle of Ulam's who "just had to go to Monte Carlo" to gamble. The first public description is:

Metropolis, N. & Ulam, S. (1949). The Monte Carlo Method. J. Amer. Statist. Assoc. 44(247):335–341.

Later Los Alamos work ran on MANIAC I (operational March 1952), the machine behind the 1953 algorithm below.

3. The problem plain Monte Carlo can't solve — and MCMC's answer

Plain Monte Carlo needs to draw independent samples from the target distribution. But the interesting distributions in physics and statistics are the ones you cannot sample from directly: they are high-dimensional, and known only up to a normalizing constant (you can compute \(\propto \pi(x)\) but not \(\pi(x)\) itself, because the normalizer is its own intractable integral).

Markov Chain Monte Carlo (MCMC) turns the problem inside out. Instead of sampling \(\pi\) directly, construct a Markov chain whose stationary distribution is \(\pi\) (stationary distributions are defined in ch. 01, §4). Run the chain; once it has mixed, the states it visits are samples from \(\pi\), and expectations under \(\pi\) become time-averages along the trajectory (the ergodic theorem for Markov chains). The samples are correlated rather than independent, but they come from the right distribution — and you never needed the normalizer.

Detailed balance is the standard way to guarantee \(\pi\) is stationary. A chain is reversible with respect to \(\pi\) if, for all states \(x, y\),

\[\pi(x)\,P(x\to y) = \pi(y)\,P(y\to x).\]

Detailed balance implies \(\pi\) is invariant, so it is the design target of almost every classical MCMC method.

4. The algorithms

Metropolis (1953). Born, again, in statistical physics — a simulation of 224 hard-disk particles to compute an equation of state, run on MANIAC I.

Metropolis, N.; Rosenbluth, A. W.; Rosenbluth, M. N.; Teller, A. H.; Teller, E. (1953). Equation of State Calculations by Fast Computing Machines. J. Chem. Phys. 21(6):1087–1092.

The rule: from state \(x\), propose a nearby move to \(x'\). With \(\pi \propto e^{-E/kT}\) (the Boltzmann distribution), let \(\Delta E = E(x') - E(x)\) and accept \(x'\) with probability

\[\min\!\left(1,\; e^{-\Delta E/kT}\right) \;=\;\min\!\left(1,\; \frac{\pi(x')}{\pi(x)}\right).\]

Downhill moves are always taken; uphill moves are taken sometimes — just often enough to satisfy detailed balance with respect to \(\pi\). (Historical note: Arianna W. Rosenbluth wrote the MANIAC implementation and Marshall Rosenbluth developed much of the method; the attribution to Metropolis is a well-documented simplification.)

Metropolis–Hastings (1970). W. K. Hastings generalized the rule to arbitrary, possibly asymmetric proposal distributions \(q(x\to x')\):

Hastings, W. K. (1970). *Monte Carlo sampling methods using Markov chains and their applications.* Biometrika 57(1):97–109.

\[\alpha(x, x') = \min\!\left(1,\; \frac{\pi(x')\,q(x'\to x)}{\pi(x)\,q(x\to x')}\right).\]

The extra ratio (the Hastings correction) restores detailed balance when the proposal is not symmetric; Metropolis is the special case \(q\) symmetric. This is the step that carried the method out of physics and into mainstream statistics.

Gibbs sampling (1984). Geman & Geman, working on Bayesian image restoration over Markov random fields, introduced updating one variable at a time from its full conditional \(\pi(x_i \mid x_{-i})\) — a special case of Metropolis–Hastings whose acceptance probability is always 1.

Geman, S. & Geman, D. (1984). *Stochastic Relaxation, Gibbs Distributions, and the Bayesian Restoration of Images.* IEEE Trans. PAMI 6(6):721–741.

Popularized for general Bayesian work by Gelfand & Smith (1990), it triggered the 1990s explosion of MCMC across statistics.

5. Where MCMC shows up now

  • Statistical physics — the founding domain: lattice models, most famously the

    Ising model (spins on a lattice, sampled from the Boltzmann distribution to find magnetization and phase transitions), lattice QCD, fluid equations of state.

  • Bayesian inference — the workhorse for drawing from a posterior

    \(p(\theta\mid \text{data}) \propto p(\text{data}\mid\theta)\,p(\theta)\) when the evidence is intractable. This is exactly the "known up to a constant" case MCMC was built for, and the reason applied Bayesian statistics was transformed after 1990.

  • Computational biology — Bayesian phylogenetics (MrBayes, BEAST), population

    genetics, sequence alignment.

  • Finance — Monte Carlo pricing of path-dependent derivatives, value-at-risk,

    stochastic-volatility estimation.

  • Modern probabilistic ML — the general engine for sampling unnormalized

    high-dimensional distributions: topic models (LDA), energy-based and generative models, and the gradient-aware descendants that scale it — Hamiltonian Monte Carlo, the No-U-Turn Sampler (NUTS, behind Stan/PyMC), and Langevin dynamics / SGLD.

Why it matters to the Stack

The Stack's data and AI layers stand on exactly the object MCMC samples: a distribution known only up to a constant, over a space too large to enumerate. Any Bayesian estimate, any probabilistic ranking, any generative sampler in a Koder component is, underneath, a Markov chain engineered to have the answer as its stationary distribution — the same construction PageRank uses to rank the web (chapter 03) and the same one language models specialize to sample text (AI Compendium, ch. 35).

Sources

  • Metropolis, N. & Ulam, S. (1949). The Monte Carlo Method. J. Amer. Statist. Assoc. 44(247):335–341.
  • Metropolis, Rosenbluth, Rosenbluth, Teller & Teller (1953). Equation of State Calculations by Fast Computing Machines. J. Chem. Phys. 21(6):1087–1092.
  • Hastings, W. K. (1970). Monte Carlo sampling methods using Markov chains and their applications. Biometrika 57(1):97–109.
  • Geman, S. & Geman, D. (1984). Stochastic Relaxation, Gibbs Distributions, and the Bayesian Restoration of Images. IEEE Trans. PAMI 6(6):721–741.
  • Eckhardt, R. (1987). Stan Ulam, John von Neumann, and the Monte Carlo Method. Los Alamos Science, Special Issue 15:131–141.
  • Gelfand, A. E. & Smith, A. F. M. (1990). Sampling-Based Approaches to Calculating Marginal Densities. J. Amer. Statist. Assoc. 85(410):398–409.