The random walk
that lands in the right places

Chapter 1: Markov chains — a memoryless walk that forgets where it started

Forget sampling for a moment. We need to understand the engine first. A Markov chain is a process that hops between states where the next state depends only on the current one — never on the path that led there. That "no memory of history" rule is the Markov property.

Take a tiny weather chain with two states. Sunny stays Sunny 70% of days (so flips to Rainy 30%); Rainy stays Rainy 60% (flips to Sunny 40%). Collect those numbers into a transition matrix P, where row from Sunny is [0.7, 0.3] and row from Rainy is [0.4, 0.6].

Start certain it is Sunny today: the distribution over states is p₀ = [1, 0]. To advance one day, multiply by P: p₁ = p₀P. Multiply again for two days, and so on. The question is: what is the long-run climate — the fraction of sunny days a year from now?

Something beautiful happens. From any starting belief — certain Sunny, certain Rainy, or 50/50 — the distribution converges to the same long-run mix. The start is forgotten. That self-forgetting limit is the chain's stationary distribution π, the one that satisfies πP = π (multiplying by P leaves it unchanged).

This convergence requires two mild conditions. The chain must be irreducible (you can eventually get from any state to any other) and aperiodic (it doesn't cycle with a fixed rhythm). Our weather chain is both, so a unique π exists and every start converges to it. (Continuous-time versions of the same idea live in CTMC, and chains with hidden states power HMMs.)

This is the lever MCMC will pull in reverse. Here you are given P and you find π. MCMC will choose the π it wants (your posterior) and engineer a P that converges to it. But first, see the forgetting happen.

Worked example: converging by hand

Transition matrix P = [[0.7, 0.3], [0.4, 0.6]], start p₀ = [1, 0]. One step:

How fast does it converge? The eigenvalues of P are always 1 (the stationary direction) and, for a 2×2 chain, λ₂ = (trace − 1) = (0.7 + 0.6) − 1 = 0.3. The distance to stationarity shrinks by a factor λ₂ = 0.3 every step — geometric convergence. Small λ₂ means fast forgetting; λ₂ near 1 means a chain that crawls (the seed of the mixing problem in Chapter 7).

Solve πP = π with π summing to 1. Writing π = [π_S, 1−π_S], balance at state Sunny gives π_S = 0.7π_S + 0.4(1−π_S), so 0.6π_S = 0.4, π_S = 4/7. Hence:

Check the forgetting: after 5 steps from [1, 0], matrix multiplication gives [0.5725, 0.4275] — essentially converged. The remaining gap is on the order of (0.3)⁵ ≈ 0.0024 of the original deviation. The chain has all but forgotten that it began certain it was Sunny.

import numpy as np
P = np.array([[0.7, 0.3], [0.4, 0.6]])
p = np.array([1., 0.])
for _ in range(50):
    p = p @ P            # repeated multiplication -> converges to pi regardless of start
print(p)                    # [0.5714 0.4286]

Or pull the stationary distribution straight out as the eigenvector of Pᵀ with eigenvalue 1:

from numpy import linalg as la
vals, vecs = la.eig(P.T)
pi = np.real(vecs[:, np.isclose(vals, 1)]).ravel(); pi /= pi.sum()  # [0.5714 0.4286]

Chapter 2: Stationarity & detailed balance — the design knob

Here is the inversion that makes MCMC possible. In Chapter 1 you were given a transition matrix and you solved for its stationary π. MCMC runs that backwards: you are handed the π you want (your target posterior) and must design a transition rule whose stationary distribution is exactly that π.

Designing for stationarity directly is awkward — πP = π is a global condition coupling every state to every other. We want something local and easy to check. That something is detailed balance.

Detailed balance says: for every pair of states x and x', the probability flux flowing x → x' equals the flux flowing back x' → x.

Read it as a plumbing condition: at equilibrium, the amount of probability "mass" leaving x toward x' is matched, pair by pair, by mass returning. A chain that satisfies detailed balance is called reversible.

Why does this give stationarity for free? Sum the detailed-balance equation over all x'. The left side becomes π(x)∑_x'P(x→x') = π(x) (the row sums to 1). The right side becomes ∑_x'π(x')P(x'→x) — exactly the total inflow to x, i.e. (πP)(x). So π(x) = (πP)(x) for every x, which is πP = π. Stationarity, delivered by a per-pair condition you can check one edge at a time.

This is the cheat code. Instead of solving a global eigenvector problem, we will build a transition rule that satisfies detailed balance edge-by-edge — and Chapter 3 shows the rule (Metropolis) does exactly that. But one warning the sim below makes vivid: detailed balance is sufficient, not necessary.

A chain can be stationary for π without being reversible. The classic counterexample is a chain that circulates probability around a loop: net flow around the cycle, yet every state's total inflow still equals its outflow, so π is unchanged. Some modern samplers deliberately break reversibility for speed — we use detailed balance because it is trivially easy to engineer, not because non-reversible chains are forbidden.

Worked example: balance on two states

Target π = [0.8, 0.2]. Use a symmetric proposal (always propose the other state) and the Metropolis accept rule P(x→x') = min(1, π(x')/π(x)). Compute the two transition probabilities.

From state 0 to 1: accept = min(1, π₁/π₀) = min(1, 0.2/0.8) = min(1, 0.25) = 0.25, so P₀₁ = 0.25.

From state 1 to 0: accept = min(1, π₀/π₁) = min(1, 0.8/0.2) = min(1, 4) = 1.0, so P₁₀ = 1.

Now the fluxes. Forward: π₀P₀₁ = 0.8 × 0.25 = 0.2. Backward: π₁P₁₀ = 0.2 × 1 = 0.2. They match, so detailed balance holds and π = [0.8, 0.2] is stationary — we engineered it, edge by edge.

Contrast the cyclic chain P = [[0.5, 0.5, 0], [0, 0.5, 0.5], [0.5, 0, 0.5]]. Its stationary distribution is uniform π = [1/3, 1/3, 1/3] (by symmetry). But check the 0↔1 edge: forward flux π₀P₀₁ = (1/3)(0.5) = 0.1667, while the backward flux π₁P₁₀ = (1/3)(0) = 0. They do not match. Detailed balance fails, yet π is still stationary — proof that balance is sufficient but not necessary.

import numpy as np
def db_holds(pi, P):
    flux = pi[:, None] * P                 # pi_i * P_ij
    return np.allclose(flux, flux.T)        # pi_i P_ij == pi_j P_ji ?

pi = np.array([0.8, 0.2])
P  = np.array([[0.75, 0.25], [1.0, 0.0]])  # Metropolis-built
print(db_holds(pi, P))                       # True

The library way is simply to construct P from π via Metropolis, which makes detailed balance automatic:

# P_ij = q_ij * min(1, pi_j / pi_i)  -> detailed balance holds by construction

Chapter 3: Metropolis — propose, compare heights, accept

Now we cash detailed balance into an actual algorithm. The Metropolis algorithm (1953) does just three things in a loop. Stand on the unnormalized target surface at point x. Propose a small symmetric nudge to a candidate x'. Then accept it with probability min(1, p(x')/p(x)) — otherwise stay put.

Unpack that acceptance rule. If the candidate is higher (more probable), the ratio p(x')/p(x) exceeds 1, min caps it at 1, and you always move. If it is lower, you move only sometimes — with probability equal to exactly how much lower it is. Always uphill; sometimes downhill.

That single asymmetry is the whole magic. Pure hill-climbing would get stuck at the mode and never sample the tails. Pure random walking would wander uniformly and ignore the shape. Metropolis threads the needle: the downhill-sometimes rule makes the walker linger in high-probability regions in the right proportion, while still visiting the tails. Its long-run visit frequency equals p.

Why does it equal p? Because the rule was reverse-engineered from Chapter 2. With a symmetric proposal, plug P(x→x') = q·min(1, p(x')/p(x)) into detailed balance and the two sides match for every pair — so p is the stationary distribution, by construction. (Chapter 4 handles the asymmetric case.)

One subtlety that trips up every beginner: a rejected proposal is not skipped. When you reject, you record the current point again as the next sample. Those repeated stays are not a bug — they are precisely how the walk spends extra time on the tall parts of the distribution. Forget to re-count them and your histogram comes out wrong.

In the sim below we target the standard normal (a stand-in posterior you can also read as a 1-D bell). Watch a translucent proposal ghost appear one symmetric step away, the height bar compare the two points, and the move accept or reject — with the histogram filling toward the true bell.

Worked example: two proposals from x = 0.5

Target: standard normal, unnormalized p(x) = exp(−x²/2). Symmetric proposal. We are at x = 0.5 and consider two candidates: x' = 1.5 (downhill, farther from the mode 0) and x' = 0.2 (uphill, toward the mode).

The acceptance ratio is p(x')/p(x) = exp(−(x'² − x²)/2). The unknown normalizer 1/√(2π) cancels — same Z-cancellation as Chapter 0.

Downhill, x' = 1.5: exponent = −(1.5² − 0.5²)/2 = −(2.25 − 0.25)/2 = −1. Ratio = exp(−1) = 0.3679. So we accept this downhill move with probability 0.3679 — taken about a third of the time.

Uphill, x' = 0.2: exponent = −(0.2² − 0.5²)/2 = −(0.04 − 0.25)/2 = +0.105. Raw ratio = exp(0.105) = 1.1107, which exceeds 1, so min(1, 1.1107) = 1.0 — always accepted.

Concretely: stepping toward the peak is automatic; stepping into the tail happens, but only with probability set by how far down the curve drops. Over millions of steps those odds reproduce the bell exactly.

import numpy as np
x, step, chain = 0.0, 1.0, []
for _ in range(100000):
    xp = x + np.random.normal(0, step)        # symmetric proposal
    a  = np.exp(-(xp**2 - x**2) / 2)        # p(xp)/p(x), Z cancels
    if np.random.rand() < a:
        x = xp                                  # accept
    chain.append(x)                              # reject re-appends current x!

The library does the same loop for you, with autotuning:

import pymc as pm
# with pm.Model(): ...; pm.sample(step=pm.Metropolis())

Chapter 4: Hastings — correcting a biased proposal

Symmetric proposals waste effort. Near a hard boundary — a rate that must stay positive, a variance that must stay above zero — a symmetric Gaussian keeps proposing illegal negative values that you must reject. You would rather use a proposal that leans toward the legal, high-probability region.

But a lopsided proposal introduces a bias: it visits some candidates more often than their mirror moves, which would distort the stationary distribution. Unless you pay a precise correction. Hastings (1970) generalized Metropolis to handle exactly this case.

Recall detailed balance: π(x)P(x→x') = π(x')P(x'→x). The transition P factors into "propose then accept": P(x→x') = q(x'|x)·A(x,x'), where q is the proposal density and A the acceptance probability. Plug that in and solve for A, and the proposal densities no longer cancel — they leave behind a ratio.

The result is the Metropolis-Hastings acceptance rule:

The first factor is the familiar target ratio. The second, the proposal ratio q(x|x')/q(x'|x), is the correction. It asks: how easy was it to propose this move, versus to propose its reverse? Moves that are "easy to make but hard to undo" get their acceptance dialed down; "hard to make but easy to undo" get it boosted. That bookkeeping restores detailed balance.

The beautiful special case: when the proposal is symmetric, q(x'|x) = q(x|x'), the proposal ratio is exactly 1 and Hastings collapses back to Metropolis. Metropolis is just Metropolis-Hastings with a symmetric proposal — you have been using the general rule all along.

In the sim, the target is a skewed positive Gamma(2,1) and the proposal leans toward the origin (asymmetric). The accept decision now shows two bars — target ratio and proposal ratio — whose product is the acceptance. Flip off the correction and watch the histogram drift to the wrong distribution.

Worked example: x = 2 proposing x' = 3

Target Gamma(2,1): unnormalized p(x) = x·e^−x for x > 0. Asymmetric proposal: given x, draw x' from an Exponential with mean x, so q(x'|x) = (1/x)e^−x'/x. We are at x = 2, proposing x' = 3.

Target ratio: p(3)/p(2) = (3e⁻³)/(2e⁻²) = 1.5 × e⁻¹ = 1.5 × 0.3679 = 0.5518. By itself this would reject the move most of the time.

Proposal ratio: q(2|3)/q(3|2). Forward q(3|2) = (1/2)e^−3/2 = 0.5 × 0.2231 = 0.11157. Reverse q(2|3) = (1/3)e^−2/3 = 0.3333 × 0.5134 = 0.17113. Ratio = 0.17113 / 0.11157 = 1.5340.

Multiply: A = min(1, 0.5518 × 1.5340) = min(1, 0.8465) = 0.8465. The proposal correction recognized that this move was "hard to propose in reverse" and boosted acceptance from a naive 0.55 up to 0.85 — exactly what keeps the chain unbiased.

import numpy as np
def p(x):  return x * np.exp(-x)            # Gamma(2,1) unnormalized
xp  = np.random.exponential(scale=x)       # asymmetric proposal, mean = x
q_fwd = (1/x) * np.exp(-xp/x)            # q(xp|x)
q_bwd = (1/xp)* np.exp(-x/xp)            # q(x|xp)
a = (p(xp)/p(x)) * (q_bwd/q_fwd)        # MH ratio = target * proposal
if np.random.rand() < a:  x = xp

Modern probabilistic programming libraries handle the proposal correction (and constrained-variable Jacobians) automatically:

# NumPyro / Stan apply the proposal + change-of-variables correction under the hood

Chapter 5: Gibbs — sample one coordinate at a time, accept always

Metropolis with one big multivariate proposal struggles as dimension grows: a single bad coordinate can sink the whole move, so acceptance plummets. But there is often a way out. You may not be able to sample the joint p(x₁, ..., x_d), yet you can sample each single coordinate given all the others.

Gibbs sampling exploits exactly this. It walks the coordinate axes one at a time. To update coordinate i, draw a fresh value from its full conditional p(x_i | x_−i) — the distribution of x_i with every other coordinate held fixed. Cycle through all coordinates, repeat.

The remarkable part: Gibbs never rejects. Its acceptance probability is identically 1. Why? Because Gibbs is a special case of Metropolis-Hastings whose proposal is the full conditional. Plug q(x_i'|x) = p(x_i'|x_−i) into the MH ratio and every factor cancels — target ratio against proposal ratio — leaving A = 1.

Intuitively, you are drawing each coordinate already "from the right distribution," so there is nothing to correct. No proposal width to tune, no rejections to waste. Gibbs is the workhorse behind classic engines like BUGS and JAGS, and it shines when the conditionals have a clean closed form (conjugate models, Gaussian fields, the Ising model of Chapter 6).

The catch is that zero rejections does not mean fast exploration. When coordinates are strongly correlated, the axis-aligned conditional moves get tiny, and the walker is forced to inch up a thin diagonal valley one short L-step at a time. The sim makes this vivid: crank the correlation toward 1 and watch the L-shaped Gibbs steps shrink to a crawl — a preview of the mixing problem.

This is the deeper lesson: "never rejects" and "mixes well" are different virtues. Gibbs has the first for free but can badly lack the second, which is why correlated posteriors often call for reparameterization or the gradient samplers of the next lesson.

Worked example: resampling x1 given x2 = 1

Target: bivariate standard normal with correlation ρ = 0.8. Its full conditionals are known in closed form: x₁ | x₂ ~ Normal(mean = ρ·x₂, variance = 1 − ρ²), and symmetrically for x₂ | x₁. Currently x₂ = 1; resample x₁.

Conditional mean: ρ·x₂ = 0.8 × 1 = 0.8. The best guess for x₁ is pulled toward 0.8 because the two coordinates are positively correlated.

Conditional variance: 1 − ρ² = 1 − 0.64 = 0.36, so the conditional standard deviation is √0.36 = 0.6. Knowing x₂ sharpened our belief about x₁ (from variance 1 down to 0.36).

So we draw x₁ ~ Normal(0.8, 0.6²). Because this draw is the exact conditional, the implicit MH acceptance is min(1, [p(new)/p(old)] × [q(old|new)/q(new|old)]) = 1 always — Gibbs accepts every proposed coordinate value.

import numpy as np
rho, x1, x2, chain = 0.8, 0.0, 0.0, []
sd = np.sqrt(1 - rho**2)                  # 0.6
for _ in range(10000):
    x1 = np.random.normal(rho*x2, sd)     # x1 | x2  (mean 0.8 when x2=1)
    x2 = np.random.normal(rho*x1, sd)     # x2 | x1
    chain.append((x1, x2))                  # no accept/reject -- always taken

The library form pulls in a Gibbs step engine (or the classic BUGS/JAGS samplers):

import pymc as pm
# pm.sample(step=pm.gibbs.Gibbs())  # Gibbs where the conditionals are tractable

Chapter 6: Showcase — the Ising spin lattice (where Metropolis was born)

The Metropolis algorithm was not invented for statistics. It was invented in 1953 at Los Alamos by Metropolis, the Rosenbluths, and the Tellers to simulate interacting particles. The cleanest instance is the Ising model: a grid of magnetic spins, each +1 or −1, that prefer to agree with their neighbors.

The physics is captured by an energy. Each pair of neighboring spins contributes −J·s_is_j: aligned neighbors (same sign) lower the energy, opposed neighbors raise it. With coupling J = 1, the system "wants" agreement. The probability of a whole configuration s is the Boltzmann distribution p(s) ∝ exp(−E(s)/T), where T is temperature.

This is a textbook target-you-can-score-but-not-normalize. The normalizer (the "partition function") sums exp(−E/T) over all 2^N configurations — astronomically intractable for any real lattice. But Metropolis needs only ratios, and a ratio of Boltzmann weights for a single spin flip depends only on the local energy change.

Flip one spin and the energy changes by dE = 2J·s_i·(sum of its four neighbors). Accept the flip with min(1, exp(−dE/T)). That is the entire algorithm: pick a spin, compute its local dE, accept or reject. Sweep this over the lattice millions of times and you are sampling configurations from the Boltzmann distribution.

The payoff is a genuine phase transition. At low temperature, energy-raising flips are almost always rejected, so spins lock into alignment — a magnet. At high temperature, even costly flips get accepted freely, thermal noise scrambles the spins, and the net magnetization collapses. The crossover happens near the critical temperature T_c ≈ 2.27 (the exact 2-D Onsager value 2/ln(1+√2)).

Drag the temperature slider below and watch order emerge and dissolve. Press "sweep T" to ramp temperature and trace the magnetization curve — you are watching a simple accept/reject rule conjure a phase transition out of nothing but local spin flips.

Worked example: two flips, two temperatures

Coupling J = 1. The energy change from flipping a spin is dE = 2J·s_i·(neighbor sum). Take a center spin s_i = +1.

Case A — all four neighbors are +1 (neighbor sum = 4): dE = 2 × 1 × (+1) × 4 = +8. At T = 2 the accept probability is exp(−8/2) = exp(−4) = 0.0183 — flipped less than 2% of the time. Aligned regions stay aligned; order is stable at low T.

Case B — three up, one down (neighbor sum = 2): dE = 2 × 1 × (+1) × 2 = +4. At T = 2.5 the accept probability is exp(−4/2.5) = exp(−1.6) = 0.2019 — flipped about a fifth of the time. These occasional thermal flips inject the disorder that erodes magnetization as T rises.

Compare the two: raising T from 2 to 2.5 (and lowering the energy cost from 8 to 4) pushed acceptance from 1.8% to 20%. That sensitivity, multiplied over a whole lattice, is what makes the transition sharp.

import numpy as np
# one Metropolis spin-flip attempt on an N x N lattice s of +/-1
nb = s[i-1,j] + s[(i+1)%N,j] + s[i,j-1] + s[i,(j+1)%N]
dE = 2 * J * s[i,j] * nb              # local energy change
if np.random.rand() < np.exp(-dE/T):
    s[i,j] *= -1                       # accept the flip

Specialized libraries ship vectorized Ising-Metropolis sweeps:

# netket and other lattice MC libraries provide optimized Metropolis sweeps over spins

No quiz this chapter — the simulation is the test. Lower T until the lattice freezes into a single color, then raise it past ~2.3 and watch order dissolve.

Chapter 7: Burn-in, mixing, and the lie of the early samples

You have to start your sampler somewhere, and that somewhere is arbitrary — maybe far out in the tail of the posterior. The first few hundred samples are then not draws from the target at all; they are the chain still hiking in from a bad starting point.

Average those early samples into your estimate and you bias it toward nowhere-land. The fix is burn-in: discard the initial stretch of the chain, keeping only samples drawn after the chain has reached its stationary regime.

How long is long enough? It depends on how fast the chain forgets its start — exactly the geometric convergence rate λ₂ from Chapter 1. The deviation from the stationary distribution shrinks by a factor λ₂ each step, so the burn-in length grows when λ₂ is near 1 (a slow, sticky chain).

After burn-in comes a second, distinct concern: mixing. Mixing is how quickly the chain moves around the target once it is there. A chain can reach stationarity (end of burn-in) yet still explore so sluggishly that consecutive samples barely differ — we tackle that in Chapter 8 with autocorrelation.

The practical tool for seeing both is the trace plot: sampled value versus iteration. Run several chains from very different starts and overlay them. While they are in burn-in the traces sit at different levels; once mixed they merge into a single fuzzy band. The simulator below lets you drag the burn-in cutoff and watch the estimate's bias appear and vanish.

A warning the trace plot also teaches: a chain that looks flat is not necessarily converged. It might be stuck (poor mixing), crawling so slowly it appears calm. Flatness alone never proves you reached the target — which is exactly why Chapter 9's multi-chain R-hat diagnostic exists.

Worked example: burn-in for the weather chain

Use the transparent 2-state chain from Chapter 1, P = [[0.7, 0.3], [0.4, 0.6]] with λ₂ = 0.3, started at [1, 0]. The deviation from the stationary π shrinks like (0.3)^k times its initial size. How many steps until we are within 1% of stationary?

We need (0.3)^k ≤ 0.01. Take logs:

Check the residual factor: (0.3)⁴ = 0.0081, comfortably under 1%. And the actual distribution after 4 steps from [1, 0] is [0.5749, 0.4251], versus the stationary [0.5714, 0.4286] — a gap of about 0.0035, well within tolerance.

So the first ~4 samples are burn-in for this fast chain. The lesson scales brutally: a slow chain with λ₂ = 0.99 would need ln(0.01)/ln(0.99) ≈ 458 steps for the same 1% — over a hundred times longer. Burn-in cost is set by how close λ₂ sits to 1.

import numpy as np
lam2 = 0.3
burn = int(np.ceil(np.log(0.01) / np.log(lam2)))   # 4  -> within 1% of stationary
samples = run_chain(N=10000)
est = samples[burn:].mean()                       # discard pre-stationary draws

In practice you eyeball the trace and let a diagnostics library confirm:

import arviz as az
az.plot_trace(idata)   # trace plots reveal burn-in and whether chains have mixed

Chapter 8: Autocorrelation and effective sample size

You collected 10,000 MCMC samples — but they are not 10,000 independent draws. Each sample is a small step from the last, so neighbors are correlated; the chain repeats itself. This is the deep difference from the plain Monte Carlo of the root lesson, where every draw was independent and the standard error was simply σ/√N.

If you naively plug N = 10,000 into σ/√N here, you will badly overstate your precision. The honest question is: how many independent samples is this correlated chain actually worth?

The answer comes from the autocorrelation function: the correlation between samples separated by a lag of 1, 2, 3, ... steps. For a well-mixed chain it crashes to zero quickly; for a sticky chain it has a long fat tail. The total area under it defines the integrated autocorrelation time τ — roughly, how many raw steps it takes to get one fresh, independent sample.

Divide your raw count by τ and you get the effective sample size:

This is the count that actually governs your Monte Carlo error. A chain of 10,000 samples with τ = 50 carries the information of only 200 independent draws.

The simulator lets you drag the proposal width to change the lag-1 correlation. Sticky chains (high ρ) grow a long autocorrelation tail and a tiny N_eff; well-mixed chains show the autocorrelation crashing to zero with N_eff near N. The two error bars — naive σ/√N and honest σ/√N_eff — are drawn side by side so you see the gap.

Worked example: an AR(1) chain with rho = 0.9

Model the chain as AR(1) with lag-1 autocorrelation ρ = 0.9 — a fairly sticky Metropolis chain. For AR(1), the autocorrelation at lag k is ρ^k, and the integrated autocorrelation time has the clean closed form τ = (1+ρ)/(1−ρ). We have N = 1000 raw samples.

So each effectively-independent sample costs 19 raw draws. The effective sample size is:

The honest standard error is σ/√52.63, which is √19 = 4.36 times wider than the naive σ/√1000. Reporting the naive error would overstate your precision by a factor of about 4.4 — the difference between a credible answer and a misleading one.

import numpy as np
rho = np.corrcoef(chain[:-1], chain[1:])[0,1]   # lag-1 autocorrelation
tau = (1 + rho) / (1 - rho)                  # AR(1) integrated autocorr time = 19
n_eff = len(chain) / tau                         # 52.6 effective draws from 1000

The library computes τ from the full autocorrelation function (more robust than the AR(1) shortcut):

import arviz as az
az.ess(idata)   # effective sample size from the autocorrelation

Chapter 9: R-hat convergence and tuning the step size

How do you know your chain converged, rather than merely looks calm? The most trusted answer runs several chains from different starting points and compares them. If they have all reached the same target, the spread between chains should match the spread within each chain.

That comparison is the Gelman-Rubin statistic R-hat. It pools a between-chain variance B (how far apart the chain means are) and a within-chain variance W (how much each chain wiggles internally) into a ratio. When the chains have mixed, R-hat → 1. When between-chain spread still dominates, R-hat sits well above 1 and you are not done.

Concretely, R-hat = √(V̂/W), where V̂ is a pooled variance estimate that blends W and B. The rule of thumb every probabilistic programming language prints: keep sampling until R-hat < 1.01 for every parameter.

The lever that controls all of this is the proposal step size. Too small, and almost every proposal is accepted but each move is microscopic — the chain crawls and mixes slowly. Too large, and bold proposals overshoot into low-probability regions and are mostly rejected — the walker gets stuck. There is a sweet spot in between.

For high-dimensional random-walk Metropolis, theory pins that sweet spot at an acceptance rate of about 0.234 — you should be rejecting roughly three out of four bold proposals. A 98% acceptance rate is a red flag, not a triumph: it almost always means your steps are far too timid.

The sim shows four chains starting spread out and slowly merging, with a live R-hat gauge ticking from large toward 1.0, alongside a step-size dial revealing the three regimes (too small, just right, too big) through their acceptance rates.

Worked example: R-hat by hand on two tiny chains

Two chains (m = 2), each of length n = 5. Chain A = [1, 2, 3, 4, 5], Chain B = [3, 4, 5, 6, 7]. Tiny and deterministic so every term is checkable.

Means: Ā = 3, B̄ = 5, grand mean = 4.

Between-chain variance: B = n/(m−1) · ∑(chain mean − grand)² = 5/1 · [(3−4)² + (5−4)²] = 5 × (1 + 1) = 10.

Within-chain variance: each chain's sample variance (ddof = 1) is 2.5 (the values 1..5 and 3..7 are just shifts), so W = (2.5 + 2.5)/2 = 2.5.

Pooled estimate: V̂ = (n−1)/n · W + B/n = (4/5)(2.5) + 10/5 = 2 + 2 = 4.

R-hat = 1.2649 is well above 1, so these chains have NOT converged — their means (3 and 5) differ too much relative to their internal spread. Keep sampling, or fix the step size, until R-hat drops below 1.01.

import numpy as np
chains = [np.array([1,2,3,4,5]), np.array([3,4,5,6,7])]
m, n = len(chains), len(chains[0])
means = np.array([c.mean() for c in chains]); grand = means.mean()
B = n/(m-1) * np.sum((means - grand)**2)          # 10
W = np.mean([np.var(c, ddof=1) for c in chains])  # 2.5
varhat = (n-1)/n * W + B/n                      # 4.0
rhat = np.sqrt(varhat / W)                       # 1.2649

In practice ArviZ prints R-hat and ESS for every parameter at once:

import arviz as az
az.summary(idata)   # r_hat and ess columns for each parameter

Chapter 10: Showcase — MCMC on combinatorial spaces (coloring & trees)

Everything so far lived on the real line, where you can take small numeric steps. But the targets that made MCMC indispensable are often discrete and non-Euclidean: which proper coloring of a graph? which evolutionary tree relates these species? You cannot take a gradient of a tree, and you cannot step a tiny ε "between" two trees.

MCMC does not care. It needs only two things: a way to propose a neighbor structure, and a way to evaluate the target ratio between the current and proposed structure. On a discrete space, the "neighbors" are defined by moves — recolor a node, swap two tree branches — and you walk a graph whose vertices are entire structures.

Graph coloring (panel A) is a clean example. Suppose you want proper 3-colorings (no two adjacent nodes share a color). Single-site Gibbs picks a node, looks at the colors its neighbors are not using, and recolors it uniformly among the allowed ones. Conflicts (same-color edges) glow red and dissolve as the chain finds valid colorings — with no gradient, no coordinates, just legal-move bookkeeping.

Bayesian phylogenetics (panel B) is where this matters most. Given DNA from several species, you want the posterior over evolutionary trees. The catch: the number of distinct unrooted binary trees on n taxa is the double factorial (2n−5)!! — it explodes super-exponentially. For just 10 species there are over two million trees; enumeration is hopeless.

So tools like MrBayes, RevBayes, and BEAST run Metropolis-Hastings over tree space. A move (a "nearest-neighbor-interchange") proposes a local rearrangement of branches; the MH ratio on the tree's likelihood accepts or rejects it. The chain concentrates on high-likelihood trees, sampling the posterior over the tree of life without ever listing all the trees.

Play with both panels below. Step the coloring until conflicts vanish; propose tree moves and watch the MH acceptance and the staggering count of possible trees the chain is navigating.

Worked example: two deterministic checks

(1) Single-site Gibbs for 3-coloring. A node has two neighbors already colored Red and Green; we use 3 colors {R, G, B}. The allowed colors for this node are those not used by neighbors: only Blue. So 1 allowed color out of 3. Single-site Gibbs draws uniformly among allowed colors, so P(Blue) = 1/1 = 1.0 — it is forced. (With a soft penalty instead of a hard constraint, it would be near-1 but not exactly 1.)

(2) The size of tree space. The number of distinct unrooted binary trees on n taxa is (2n−5)!!. For n = 5: (2×5−5)!! = 5!! = 5 × 3 × 1 = 15 trees. For n = 10: (2×10−5)!! = 15!! = 15·13·11·9·7·5·3·1 = 2,027,025 trees.

Two million trees for ten species — and real studies have hundreds of taxa, pushing the count past the number of atoms in the universe. Enumeration is impossible; MCMC's move-by-move walk is the only route.

import random
# single-site Gibbs recolor: draw uniformly among colors not used by neighbors
allowed = [c for c in colors if c not in neighbor_colors(node)]
node_color[node] = random.choice(allowed)   # no gradient needed

# size of unrooted binary tree space on n taxa: (2n-5)!!
def n_trees(n):
    p = 1
    for k in range(2*n-5, 0, -2): p *= k
    return p
n_trees(10)   # 2027025

Production phylogenetics packages do this at scale:

# MrBayes / RevBayes / BEAST run Metropolis-Hastings over tree space for Bayesian phylogenetics

No quiz this chapter — the simulation is the test. Step the coloring until zero conflicts remain, then propose tree moves and watch the search space size you are navigating.

Chapter 11: Connections — where MCMC goes next

You now own the core. Build a chain whose stationary distribution is your target, average along it (after burn-in), and diagnose convergence with R-hat and effective sample size. That toolkit alone powers an enormous slice of modern Bayesian statistics.

But random-walk Metropolis has one structural weakness: it proposes blindly. In high dimensions, a blind step is almost always a bad step, so the optimal step size shrinks and mixing slows down. The frontier asks: what if the proposal used the gradient of the log-target to propose intelligently, like a ball rolling on the energy surface?

That is exactly the next lesson, Hamiltonian Monte Carlo. HMC and its auto-tuned cousin NUTS, plus Langevin dynamics, use the slope of the log-density to glide along the distribution's geometry instead of stumbling across it. They are what Stan, PyMC, and NumPyro actually run under the hood — and the reason Bayesian deep learning at scale is even possible.

The sim below races a random-walk walker against a gradient-guided one on the same correlated target. Raise the dimension and watch the gap widen — the quantitative case for why the next lesson exists.

Beyond gradients, MCMC is a whole family. Simulated annealing, parallel tempering, and sequential Monte Carlo all build chains for hard targets; continuous-time chains and hidden Markov models share the Markov machinery you learned in Chapter 1; and MCMC inference over latent variables is the engine inside many decision and control methods. The principle never changes: engineer a chain whose stationary distribution is your target.

And it is genuinely current. The NUTS engines in Stan / PyMC / NumPyro power today's Bayesian workflows; posterior sampling drives uncertainty quantification in Bayesian deep learning; and MrBayes / RevBayes / BEAST reconstruct the tree of life from genomic data — all built on the accept/reject rule you just learned by hand.

Worked example: why dimension forces the issue

For an isotropic Gaussian target in d dimensions, the optimal random-walk step scales like 1/√d, so the per-step squared distance — and hence the effective samples per step — degrades roughly as 1/d. Compare the relative cost at d = 1, 10, 100.

Take ESS-per-step proportional to 1/d. Relative to d = 1: at d = 10 efficiency is 1/10 = 0.1; at d = 100 it is 1/100 = 0.01.

So a 100-dimensional posterior needs roughly 100× more random-walk steps per independent sample than a 1-D one. Well-tuned HMC degrades far more gently (closer to d^1/4 cost), which is precisely why gradient samplers took over. (These are illustrative scalings, not exact constants.)

# random walk:           xp = x + step * randn(d)
# gradient-guided (Langevin preview): xp = x + 0.5*eps**2 * grad_logp(x) + eps * randn(d)

import numpyro
numpyro.infer.MCMC(numpyro.infer.NUTS(model)).run(key)   # gradient MCMC with autodiff