A problem-driven introduction to Monte Carlo, MCMC, and related sampling methods.
Sampling algorithms are usually presented as a sequence of methods. This tour presents them as a sequence of problems. Each chapter identifies a limitation of the preceding method, demonstrates it in the interactive gallery, quantifies how it scales with dimension using experiments that run directly in the browser, and then introduces the method designed to address it. The aim is to make clear not only how each algorithm works, but why it is needed.
Chapter 1 · Numerical integration in high dimensions
The integration problem
Bayesian inference — and much of applied probability — reduces to computing expectations,
Eπ[f] = ∫ f(x) π(x) dx.
The oldest numerical answer is deterministic quadrature: lay a grid of m points along each axis and sum. In one dimension this works well — the midpoint rule's error falls like m−2. But a grid over d dimensions requires n = md evaluations, so at a fixed budget n each axis receives only n1/d points and the error scales as
error ∝ n−2/d.
By d = 6, a budget of several thousand evaluations provides only three points per axis. This is the curse of dimensionality.
Monte Carlo takes a different approach: rather than covering the space with a grid, it averages f over random draws x⁽¹⁾,…,x⁽ⁿ⁾ ~ π,
The error of the average is governed by the central limit theorem, and the dimension d appears nowhere in the rate.
The rates
Quadrature (midpoint, product rule): error ∝ n−2/d — the exponent deteriorates with dimension.
Monte Carlo: error = sd(f(X))·n−1/2 — the same rate in every dimension; d enters only through the constant.
Refs: Davis & Rabinowitz, Methods of Numerical Integration; Robert & Casella, Monte Carlo Statistical Methods, ch. 3.
Interlude — the distance from the origin of a d-dimensional Gaussian
The distance from the origin of Gaussian draws, resampled continuously. The density is always maximised at the origin, yet as d increases the samples concentrate in a thin shell near radius √d (dashed line). In high dimension the point of highest density and the region containing most of the probability mass are far apart; this fact underlies the failures examined in the following chapters.
Experiment — estimating E[‖x‖²/d] = 1 under a Gaussian at a fixed evaluation budget
Computed live in the browser. Quadrature is far more accurate while the grid is dense, but its error grows rapidly once n1/d falls to a few points per axis, while the Monte Carlo error is nearly independent of d. Increasing the budget sixteen-fold moves the crossover point by roughly one dimension. Log scale: each gridline is a factor of ten.
Monte Carlo's guarantee presumes the ability to draw samples from π; producing those samples is itself a problem. → Chapter 2.
Chapter 2 · Problem: Monte Carlo requires samples from π
Generating exact samples: rejection sampling
For a handful of textbook distributions, samples can be generated by inverting the distribution function. A Bayesian posterior, however, is typically known only up to a normalising constant, π(x) ∝ γ(x), and no direct formula produces draws from it.
Rejection sampling is the classical answer. Draw x ~ q from a proposal that can be sampled directly, and accept it with probability
α(x) = γ(x) / (M q(x)), where M ≥ supx γ(x)/q(x).
The envelope constant M guarantees that α(x) ≤ 1, and the accepted points are exact, independent samples from π. For normalised γ the overall acceptance rate is 1/M, so the entire cost of the method is concentrated in M.
This failure is not an artefact of poor tuning; it is the curse of dimensionality in another form. The acceptance rate is 1/M, and M is the worst-case ratio of target to proposal over all of space. In high dimension, mild mismatch compounds coordinate by coordinate.
The rate
Target N(0, Id), proposal N(0, σ²Id) with σ > 1: the tightest envelope is M = σd, so the acceptance rate is σ−d — exponentially small in d. At σ = 1.5, dimension 16 accepts roughly one draw in 660.
Refs: von Neumann (1951); MacKay, ITILA, ch. 29.
Experiment — acceptance rate vs dimension, 40,000 proposals per point
Dots: measured acceptance. Dashed: the theoretical rate σ−d. Reducing σ towards 1 slows the decay, but the decay disappears only at σ = 1 — that is, when the proposal coincides with the target.
Each rejected draw is discarded computation. Importance sampling retains every draw and accounts for the mismatch through weights. → Chapter 3.
Chapter 3 · Problem: rejection sampling discards most of its draws
Weighting instead of discarding: importance sampling
Importance sampling retains every proposal draw and corrects for the mismatch with a weight. Expectations become weighted averages,
Ê[f] = ∑i wi f(xi) / ∑i wi, wi = γ(xi)/q(xi).
Nothing is discarded, no envelope is required, and the estimator is consistent for any q whose support covers the target. The cost is subtler: the draws are no longer equally informative. A small number of points with very large weights can dominate the sample — weight degeneracy — and the appropriate measure of sample quality is the weighted effective sample size,
ESS = (∑i wi)² / ∑i wi²,
which equals n when the weights are uniform and approaches 1 when a single weight dominates.
The curse of dimensionality reappears through the variance of the weights. Even a well-placed, safely over-dispersed proposal loses effective samples at a geometric rate as the dimension grows — mismatch multiplies across coordinates exactly as rejection's envelope did.
The rate
Target N(0, Id), proposal N(0, σ²Id): the effective sample fraction is ESS/n = (σ⁴/(2σ²−1))−d/2 — exponential decay in d for any σ ≠ 1. At σ = 1.5, dimension 32 keeps two effective samples per thousand.
Refs: Kish (1965); Agapiou, Papaspiliopoulos, Sanz-Alonso & Stuart (2017), Statist. Sci.; Vehtari et al. (2024), JMLR (PSIS).
Experiment — weighted ESS fraction vs dimension, 20,000 draws per point
Dots: measured Kish ESS ÷ n. Dashed: the closed-form rate. Changing σ alters the constant but not the exponential form of the decay: no fixed proposal remains effective in high dimension.
The common cause of these failures is the reliance on a single, fixed, global proposal. The alternative is to propose locally from the current state, at the cost of producing correlated samples: Markov chain Monte Carlo. → Chapter 4.
Chapter 4 · Problem: any fixed global proposal inherits the curse
From independent draws to Markov chains: Metropolis–Hastings
Rejection and importance sampling failed for one shared reason: a single global proposal q must match π everywhere at once, and per-coordinate mismatch multiplies into an exponential penalty. Metropolis–Hastings removes that requirement by proposing locally and constructing a Markov chain whose stationary distribution is exactly π.
One step of the random-walk algorithm, from the current state x:
propose x′ = x + σξ, with ξ ~ N(0, I);
accept x′ with probability α = min(1, γ(x′)/γ(x));
otherwise the next state remains x, and x is recorded again.
Only density ratios appear, so the unnormalised γ suffices; no envelope is required, and there is no global proposal to tune. The cost is correlation: successive states are dependent, so n steps of the chain are worth fewer than n independent samples. The deficit is quantified by the effective sample size,
ESS = n / (1 + 2 ∑k≥1 ρk),
where ρk is the lag-k autocorrelation of the chain — the quantity reported throughout the gallery.
On N(0, Id) with the optimally scaled proposal σ = 2.38/√d, the acceptance rate tends to the well-known value 0.234 and the efficiency per step decays like 1/d — polynomial in d, where rejection and importance sampling were exponential. This exchange of exponential for polynomial cost is the reason MCMC is the standard tool in high dimensions.
Refs: Metropolis et al. (1953); Hastings (1970); Roberts, Gelman & Gilks (1997), Ann. Appl. Probab.
Experiment — ESS fraction vs dimension: importance sampling against random walk Metropolis
Teal dots: random walk Metropolis with optimal scaling, measured. Grey dashed: a 1/d reference. Orange dashed: chapter 3's importance-sampling rate on the same axes. Both decay, but one polynomially and the other exponentially.
A 1/d decay is a substantial improvement, but the random walk ignores information that is usually available at little extra cost: the gradient of the log-density. → Chapter 5.
Chapter 5 · Problem: random-walk exploration is diffusive and slow
Using gradient information: MALA and Hamiltonian Monte Carlo
A random walk requires roughly (L/σ)² steps to travel a distance L. The gradient of log π indicates the direction of increasing density and is usually cheap to evaluate relative to the density itself; the methods of this chapter put it to work.
MALA adds the gradient as a drift on each proposal,
x′ = x + (σ²/2) ∇log π(x) + σξ, ξ ~ N(0, I),
followed by a Metropolis–Hastings correction with the (asymmetric) proposal density. Hamiltonian Monte Carlo goes further: the state is augmented with a momentum p ~ N(0, I) and the pair evolves under Hamiltonian dynamics for the energy
H(x, p) = −log π(x) + ½‖p‖²,
simulated with L leapfrog steps of size ε. Successive positions lie on a single coherent trajectory rather than a diffusive path. The benefit appears in the step-size exponents: to hold the acceptance rate fixed as d grows, the random walk must scale σ ∝ d−1/2, MALA only σ ∝ d−1/6, and HMC ε ∝ d−1/4, while producing nearly independent samples along each trajectory.
Evaluations needed per effective sample: random walk ∝ d², MALA ∝ d4/3, HMC ∝ d5/4 — equivalently, efficiency-per-evaluation slopes of −1, −1/3 and −1/4 on a log-log plot. The lab below measures all three.
Refs: Roberts & Rosenthal (1998, 2001); Neal (2011), Handbook of MCMC; Beskos et al. (2013), Bernoulli.
Experiment — ESS per 1000 target/gradient evaluations, each sampler optimally scaled
Measured at each dimension; the steps in the HMC curve arise because its leapfrog count is an integer. The log-log slopes are approximately −1, −1/3 and −1/4, in agreement with the theory.
This efficiency depends on two tuning parameters — the step size ε and the path length L — and it is sensitive to both. → Chapter 6.
Chapter 6 · Problem: HMC is sensitive to its tuning parameters
Tuning and robustness
As ε increases, HMC's efficiency improves until the leapfrog integrator crosses its stability threshold, at which point the acceptance rate collapses. The path length is equally delicate: efficiency is cyclic in the integration time T = Lε, since trajectories that double back waste their gradient evaluations.
The gallery contains three responses to this problem:
NUTS tunes both parameters automatically — doubling the trajectory until it U-turns, and adapting ε by dual averaging. It is the engine of Stan.
AAPS is constructed so that its tuning parameter matters far less: it proposes from a whole window of apogee-to-apogee segments, and its acceptance rule uses weighted sums rather than a single endpoint energy.
In the MALA regime, the Barker proposal trades some peak efficiency for robustness to step-size misspecification.
The leapfrog is stable only while ε < 2/√λmax (the target's sharpest curvature); beyond it the energy error grows without bound and everything is rejected. NUTS locates the efficient region automatically; AAPS reduces the sensitivity — its weighted-sum acceptance degrades gradually as ε grows, and its efficiency is nearly flat in K.
Experiment — efficiency vs step size on the banana: the effect of L and K
ESS per 1000 gradient evaluations, averaged over five seeds. HMC rises to a plateau and then falls sharply once ε exceeds its stability threshold (≈ 0.45 here); AAPS continues to improve at ε = 0.5 and degrades more gradually. Beyond ε ≈ 0.55 the integrator itself fails on this target. NUTS is not plotted, as it selects the step size automatically. Varying L reshapes the HMC curve substantially; varying K leaves the AAPS curve almost unchanged.
All of the above assumes that a single step size is appropriate throughout the space. For some targets no such step size exists. → Chapter 7.
Chapter 7 · Problem: no single step size suits the whole target
Position-dependent geometry
Neal's funnel couples a scale variable to the coordinates it governs. In the gallery's coordinates, the vertical variable y sets the scale of the horizontal one,
y ~ N(2, 3²), x | y ~ N(0, ey−2),
so the mouth of the funnel is orders of magnitude wider than the neck, and both are part of the same posterior — exactly the geometry of hierarchical models. A step size small enough for the neck wastes evaluations in the mouth; one sized for the mouth diverges in the neck. Fixed-metric HMC does not merely slow down: it under-samples the neck and produces biased estimates.
Riemannian-manifold HMC replaces the fixed mass matrix with a position-dependent metric G(x) built from local curvature, giving the non-separable Hamiltonian
H(x, p) = −log π(x) + ½ log det G(x) + ½ pᵀG(x)⁻¹p.
The sampler's notion of distance adapts to the target, so a single dimensionless step size is appropriate everywhere. The gallery draws G(x) as ellipses along the trajectory; in the neck they become strongly elongated.
Measured with the gallery's own implementations (3000 steps; the funnel's vertical marginal is exactly N(2, 3))
sampler
mean
sd
exact
2.00
3.00
HMC, tuned defaults
2.31
2.85
Riemannian HMC (SoftAbs)
1.96
3.04
The computational cost is real: each RMHMC step solves implicit equations involving Hessians. Refs: Neal (2003); Girolami & Calderhead (2011), JRSS B; Betancourt (2013), GSI.
Two further problems remain; both are active areas of research. → Chapter 8.
Chapter 8 · Two remaining problems
Multimodality and large datasets
Modes. Every sampler in chapters 4–7 is local, and a local sampler crosses a region of low probability at a rate that is exponentially small in the barrier height — it may be formally ergodic yet remain in one mode for the whole of a practical run. The remedy is to temper: work with the family of flattened distributions
πβ(x) ∝ π(x)β, 0 < β ≤ 1.
Parallel tempering runs a ladder of these copies in parallel and exchanges states between adjacent levels, so that barrier crossings made at small β propagate to the target at β = 1. Tempered SMC instead moves a weighted population along the same path from β ≈ 0 to β = 1. Both discover modes systematically rather than by chance.
Experiment — time spent in the right-hand mode of an equal 50/50 mixture, target 0.5
Running occupancy of the right-hand mode, drawn in step order; both samplers start in the left mode. At small separation the random walk crosses between modes freely. As the modes move apart its crossing rate falls exponentially — at ±4 it typically completes all six thousand steps without a single crossing — while parallel tempering continues to converge to 0.5.
Data. Every Metropolis correction in this tour evaluates π — in Bayesian practice, a product over the whole dataset — at every step. For large datasets this evaluation is the dominant cost.
The current research response is to subsample. Stochastic gradient Langevin dynamics runs the diffusion of chapter 5 on minibatch gradients with no correction step, and stochastic-gradient PDMPs drive non-reversible dynamics in the same way; in both cases the bias is controlled by the step size or the gradient refresh window. Both are in the gallery, where the gradient noise can be varied directly: see SGLD and the stochastic-gradient PDMPs, and their field notes for the details.
Summary
Eight chapters: quadrature → Monte Carlo → rejection → importance weights → Markov chains → gradients → self-tuning → geometry → populations and minibatches. Each method was introduced in response to a demonstrated limitation of the one before. The gallery contains every sampler discussed here, and the field notes provide further detail on each.
Refs: Geyer (1991); Del Moral, Doucet & Jasra (2006), JRSS B; Welling & Teh (2011), ICML; Fearnhead, Grazzi, Nemeth & Roberts (2024), arXiv:2406.19051.