Learn AI Create the Future — with Pranay
Main site

Resampling, Likelihood and Bayes

MCMC from scratch

Markov chain Monte Carlo explores a probability distribution by wandering through it — build the Metropolis algorithm in twenty lines and check it against an exact answer.

Read these first

On this page 5
  1. Why it exists
  2. How it works
  3. A real example you have seen
  4. Remember this
  5. What to learn next

One lesson, three depths. Pick the one that fits you today — you can switch any time.

Beginner — No maths. Plain English.

MCMC draws samples from a complicated distribution by taking a long random walk that lingers where probability is high.

A street-food blogger explores a new city with a rule. Each evening she considers a random nearby stall. If it is busier than tonight's, she moves there tomorrow. If quieter, she sometimes moves anyway — small risks keep the tour honest.

Follow her for a year and count the nights: her diary lists every neighbourhood in proportion to how good its food is. Nobody gave her a map. The walk became the map.

Why it exists

The conjugate shortcut gives exact posteriors — for a lucky handful of models. Step outside that handful and the belief distribution you want has no formula anyone can integrate.

Here is the strange loophole. For most models you can cheaply score any single candidate: "how well does this parameter value explain the data?" You still cannot describe the whole landscape. MCMC turns that one weak ability into full exploration. It was invented for hydrogen-bomb physics in 1953; statisticians realised in 1990 that it unlocked almost every Bayesian model, and the field exploded.

The name, unpacked: Monte Carlo = answering by random sampling. Markov chain = a walk whose next step depends only on where you stand now.

How it works

stand at a parameter value
  → propose a small random step
  → score both spots: how well does each explain the data?
      uphill (better score)?    take the step
      downhill (worse score)?   take it sometimes —
                                the worse it is, the rarer
  → repeat 20,000 times

nights spent at each value  ≈  your belief in that value

The "sometimes accept worse" rule is the soul of the method. Always climbing would strand you at one peak. Occasional descents let the walk cross valleys and give every region its fair share of visits — exactly its share, in the long run.

The early steps still remember the arbitrary starting point, so everyone discards them. That throwaway phase is called burn-in.

A real example you have seen

Election-forecast models (like FiveThirtyEight's) run MCMC to combine hundreds of polls into seat probabilities. Hurricane-track cones, drug-dose safety curves, and the black-hole image reconstruction all ran on it. Under the hood, tools named Stan and PyMC carry this walk into thousands of dimensions.

Remember this

  • MCMC needs only a scorer — no formula for the whole distribution.
  • Uphill steps always; downhill steps sometimes — that balance makes visit-counts match probability.
  • Discard the burn-in; the rest of the walk is your posterior.

What to learn next

Developer — Code and libraries.

Setup

bash
pip install numpy scipy

Outputs verified with numpy 1.26 and scipy 1.14, CPU. The chain reproduces with these versions; other versions may shift the third decimal.

Metropolis in twenty lines, with an answer key

The problem: a feature succeeded 9 times in 12 demos; we want the posterior over its success rate (flat prior). We know the exact answer from conjugacy — Beta(10, 4) — which makes this the perfect test track: build the walker, then check it against truth.

metropolis.py
import numpy as np
from scipy import stats

rng = np.random.default_rng(8)

def log_posterior(p):
    # Score for candidate rate p: 9 successes, 3 failures, flat prior.
    if not 0 < p < 1:
        return -np.inf                     # outside (0,1): impossible
    return 9 * np.log(p) + 3 * np.log(1 - p)

samples, p = [], 0.5                       # start anywhere reasonable
for _ in range(20_000):
    proposal = p + rng.normal(0, 0.1)      # a small random step
    if np.log(rng.uniform()) < log_posterior(proposal) - log_posterior(p):
        p = proposal                       # accept: always if uphill,
    samples.append(p)                      # sometimes if downhill

chain = np.array(samples[2000:])           # drop the burn-in
exact = stats.beta(10, 4)
print(f"MCMC posterior mean:  {chain.mean():.3f}")
print(f"exact answer:         {exact.mean():.3f}")
print(f"MCMC 95% interval:  ({np.quantile(chain, .025):.3f}, {np.quantile(chain, .975):.3f})")
print(f"exact 95% interval: ({exact.ppf(.025):.3f}, {exact.ppf(.975):.3f})")
Output
MCMC posterior mean:  0.715
exact answer:         0.714
MCMC 95% interval:  (0.465, 0.907)
exact 95% interval: (0.462, 0.909)

A random walk, knowing nothing but a scoring function, recovered the true posterior to about two decimal places.

The walkthrough

The acceptance line is the whole algorithm. log(uniform) < score_new - score_old accepts with probability min(1, ratio_of_posteriors). Uphill: the ratio exceeds 1, always accepted. Downhill: accepted with probability equal to the ratio — the blogger's "sometimes".

Everything happens in logs. Posterior ratios become score differences, and tiny probabilities never underflow. Same discipline as maximum likelihood.

The scorer ignores constants. The true posterior includes a normalising constant we never computed — the acceptance ratio cancels it. That cancellation is why MCMC exists: the impossible integral drops out of the algebra.

Step size 0.1 is a tuned choice. Steps of 0.001 accept nearly always but crawl; steps of 2.0 propose absurdities and stall in place. Folklore target: accept roughly a quarter to half of proposals. Track it with a counter when you experiment.

samples.append(p) runs even on rejection — staying put is a legitimate visit and the counts are wrong without it. This is the most common from-scratch implementation bug.

Common mistakes

Keeping the burn-in. The first stretch documents your starting guess, not the posterior. Drop it; for real work, run several chains from scattered starts and confirm they agree (the R-hat diagnostic).

Treating the chain as independent draws. Consecutive samples are correlated — 18,000 steps may hold the information of a few hundred independent ones (the effective sample size). Uncertainty-of-the-estimate calculations must use ESS, not chain length.

Forgetting the bounds check. Without the -np.inf guard, a proposal of 1.02 feeds log(1 - 1.02) and crashes with RuntimeWarning: invalid value encountered in log. Guard the support explicitly.

Hand-rolling Metropolis for serious multi-dimensional work. This lesson's walker is for understanding. Real problems use gradient-guided samplers (Hamiltonian Monte Carlo / NUTS) via PyMC or Stan — same idea, far fewer steps wasted.

Try it yourself

Add an acceptance counter and rerun with step sizes 0.01, 0.1, and 1.0. Record acceptance rate and the estimated mean's error for each. You have reproduced the tuning trade-off every MCMC practitioner lives with.

What to learn next

Researcher — Mathematics and papers.

The Metropolis–Hastings kernel

Target $\pi(\theta)$ known up to a constant. From state $\theta$, propose $\theta' \sim q(\theta' \mid \theta)$ and accept with probability

$$ \alpha(\theta, \theta') = \min!\left(1, \; \frac{\pi(\theta')\, q(\theta \mid \theta')}{\pi(\theta)\, q(\theta' \mid \theta)}\right) $$

Where:

  • $q$ — the proposal density (symmetric $q$ cancels, giving the Metropolis form used above).
  • On rejection the chain repeats $\theta$.

The resulting kernel satisfies detailed balance, $\pi(\theta) P(\theta \to \theta') = \pi(\theta') P(\theta' \to \theta)$, hence $\pi$ is stationary. With irreducibility and aperiodicity, the ergodic theorem gives, for integrable $f$:

$$ \frac{1}{N} \sum_{t=1}^{N} f(\theta_t) \;\xrightarrow{a.s.}\; \mathbb{E}_\pi[f] $$

A CLT holds under geometric ergodicity, with asymptotic variance $\sigma_f^2 = \operatorname{Var}_\pi(f) \cdot \tau_f$, where $\tau_f = 1 + 2\sum_{k\geq1} \rho_k$ is the integrated autocorrelation time; effective sample size is $N / \tau_f$.

Tuning theory

For random-walk Metropolis on product targets in high dimension, the optimal acceptance rate is $0.234$, with proposal scale $\propto 2.38\, \sigma / \sqrt{d}$ (Roberts, Gelman and Gilks, 1997) — mixing time grows as $O(d)$ in the step-count sense. Gradient-based samplers improve the scaling: MALA (Langevin proposals) targets $0.574$ acceptance and scales as $O(d^{1/3})$; Hamiltonian Monte Carlo simulates Hamiltonian dynamics with leapfrog integration, suppressing random-walk behaviour and scaling near $O(d^{1/4})$ (Neal, 2011, MCMC using Hamiltonian dynamics; Duane et al., 1987). NUTS (Hoffman and Gelman, 2014) removes HMC's trajectory-length tuning and powers Stan and PyMC defaults.

Other kernels

  • Gibbs sampling (Geman and Geman, 1984): cycle through full conditionals $\theta_j \sim \pi(\theta_j \mid \theta_{-j})$ — acceptance-free; a special case of MH with $\alpha \equiv 1$.
  • Data augmentation (Tanner and Wong, 1987): introduce latents to make conditionals tractable.
  • Parallel tempering: run chains at flattened powers $\pi^{1/T}$ and swap — the standard remedy for well-separated modes, which plain Metropolis crosses exponentially rarely.
  • Reversible-jump MCMC (Green, 1995) moves across models of different dimension.

Diagnostics

No finite diagnostic proves convergence; the standard battery: multiple over-dispersed chains with rank-normalised split-$\hat{R} < 1.01$ (Vehtari et al., 2021), ESS per parameter, trace plots, and divergence counts for HMC (which localise regions the sampler cannot enter — often signalling a pathological posterior geometry, cured by reparameterisation).

History and complexity

Metropolis, Rosenbluth, Rosenbluth, Teller and Teller (1953, Equation of state calculations by fast computing machines) — hard-sphere physics on the MANIAC. Hastings (1970) generalised to asymmetric proposals; Gelfand and Smith (1990, Sampling-based approaches to calculating marginal densities) detonated the Bayesian revolution. Cost per step is one posterior evaluation — for big data, subsampling variants (stochastic-gradient MCMC: Welling and Teh, 2011) trade exactness for scalability, and variational inference abandons sampling entirely for optimisation.

What to learn next