Workshop 02 — Sampling: How the Machine Works

In workshop 00 you computed a posterior the honest way: a fine grid, prior × likelihood at every point, normalize. For one unknown that works beautifully. An MMM has dozens of unknowns — and the grid trick dies instantly. This page opens the hood on what PyMC does instead: it never maps the posterior, it draws samples from it. You'll build the sampler, drive it, break it with the step-size slider, and then watch R-hat and ESS — computed live on the very chains you generated — catch every failure you just caused.

1 — Why grids die: count the cells

A small MMM has a coefficient, a carryover rate, and a saturation parameter per channel, plus trend, seasonality, and noise terms — call it 40 parameters. Lay a grid with 100 points per axis over that and you need \(100^{40} = 10^{80}\) evaluations of prior × likelihood. This explosion is the curse of dimensionality — every grid-style method becomes impossible as the number of parameters grows, because the number of grid cells multiplies with every new parameter. Suppose one evaluation takes a nanosecond — a generous one billion evaluations per second. Count it yourself:

Computed live in your browser — the same math as the notebook. Each bar is the wall-clock time for a full grid at 1 billion evaluations/second, on a log scale spanning 70+ orders of magnitude. The notebook's table pins the punchlines at 100 points/axis: 5 parameters ≈ 10 s, 10 parameters ≈ 3,170 years, 40 parameters ≈ 2.3 × 10⁵³ times the age of the universe. The grid that took microseconds in workshop 00 is physically impossible at MMM scale.

2 — Monte Carlo: answer questions by sampling

The escape hatch: instead of mapping a quantity everywhere, estimate it from random samples. This is the Monte Carlo method — replace an exact calculation with the average of many random draws — and its magical property is that its cost doesn't depend on the number of dimensions, only on how many samples you take.

The classic demo: estimate π by throwing darts uniformly at a 1×1 square. The quarter-circle inside has area \(\pi/4\), so \(\hat{\pi} = 4 \times \frac{\text{darts inside}}{\text{darts thrown}}\). No geometry, no calculus — just counting. Drag the slider and watch the estimate sharpen:

Computed live in your browser — the same math as the notebook (seeded mulberry32 instead of numpy's PCG64, so the exact digits differ run to run; the notebook's seeded 5,000-dart run landed at π̂ = 3.1728 vs truth 3.1416). Left: green darts in, red darts out. Right: the running estimate converging on π as darts accumulate.

How fast does the answer improve? Repeat the experiment many times per dart-count and measure the typical error. On a log–log plot it falls on a straight line of slope ≈ −½: Monte Carlo error shrinks like \(1/\sqrt{n}\). That's the root-n rate — to cut the error in half, you need four times the samples. Slow — but it never gets slower in higher dimensions, which is the whole trade.

Computed live in your browser — 100 repeat experiments per dart count, seeded. The notebook's heavier version (200 repeats) measured a log–log slope of −0.49 ≈ −0.5, with typical error shrinking 16× from n = 100 to n = 30,000. Your live slope is printed on the chart — it should land near −0.5 too.

3 — One catch: you can't throw darts at a posterior

The dart trick worked because we know how to draw uniform points in a square. But the thing we want to sample is the posterior — and for any real model, nobody knows how to draw from it directly. From workshop 00, Bayes' rule says

\[ p(\theta \mid \text{data}) = \frac{p(\text{data} \mid \theta)\, p(\theta)}{p(\text{data})}. \]

The numerator — likelihood × prior — is easy: pick any candidate \(\theta\), plug in, get a number. The denominator \(p(\text{data})\) is the killer: it's the sum of the numerator over all possible \(\theta\) — exactly the universe-sized grid we just gave up on. So:

• ✅ we can compute a function proportional to the posterior at any point we choose,
• ❌ we cannot compute the normalizing constant, and we cannot draw samples directly.

The solution is MCMC — Markov chain Monte Carlo: send a walker wandering through parameter space, with movement rules rigged so that, in the long run, the walker spends time in each region in exact proportion to its posterior probability. The walker's visit log is a sample from the posterior — and the rules only ever need the computable numerator. The unknowable denominator cancels.

The notebook makes this concrete with a target it knows the answer to: an email campaign got 12 signups from 100 visitors, prior Beta(2, 8), so the exact posterior is Beta(14, 96) with mean 0.1273 and sd 0.0316 — an answer key to grade the sampler against.

Computed live in your browser — the same math as the notebook. Left: what MCMC can evaluate — prior × likelihood, unnormalized (its peak height here is ~10⁻¹⁸; the scale is meaningless). Right: the normalized posterior we actually want. Same shape, unknown scale — and sampling only needs the shape.

4 — Build a sampler in ~20 lines, then drive it

The oldest MCMC recipe (Metropolis et al., 1953 — it predates the moon landing) is short enough to memorize. The walker stands at position \(x\) and repeats:

1. Propose a move: a candidate \(x' = x + \text{noise}\), e.g. a Gaussian step. The candidate-generating rule is the proposal — how the walker suggests its next location.
2. Compare heights: \(r = \tilde p(x') / \tilde p(x)\), the ratio of unnormalized posterior heights. The unknown constant divides out — this is the trick.
3. Accept or reject: uphill (\(r \ge 1\)) always moves; downhill moves with probability \(r\) — flip a biased coin. If rejected, stay put and record the current position again. The fraction of proposals taken is the acceptance rate.

Rule 3 is the genius part: always rising would find the peak and stop (that's optimization); accepting downhill moves in proportion to relative height makes time spent proportional to posterior height. The recorded list of positions is a chain — the walker's visit log, which becomes your bag of posterior samples. Here is the notebook's entire sampler:

def metropolis(log_target, start, step, n, seed):
    rng = np.random.default_rng(seed)
    chain = np.empty(n)
    x, lx, n_acc = float(start), log_target(start), 0
    for i in range(n):
        prop = x + rng.normal(0.0, step)        # 1. propose
        lp = log_target(prop)                   # 2. compare heights (in logs)
        if np.log(rng.uniform()) < lp - lx:     # 3. accept downhill with prob r
            x, lx, n_acc = prop, lp, n_acc + 1
        chain[i] = x                            # rejected -> record x again
    return chain, n_acc / n

On the Beta(14, 96) target, the notebook's seeded run of 20,000 steps took 0.04 s, accepted 44% of proposals, and recovered sample mean 0.1270 vs the exact 0.1273 and sample sd 0.0313 vs exact 0.0316 — and its 90% interval [0.080, 0.182] matched the exact [0.079, 0.183]. Twenty lines of numpy, knowing only an unnormalized log-height, reproduced the exact posterior.

🎛️ Drive the walker yourself

Below is the same algorithm, implemented in JavaScript and running in your browser, on a 2-D posterior — a correlated ridge (the long diagonal valley that real MMM posteriors are full of). Four walkers start in the four corners of the wasteland. Press ▶ Animate to watch them step: green diamonds are accepted proposals, red diamonds rejected ones (the walker stays put). Then grab the step-size slider and break it:

• Tiny steps — almost everything is accepted (nearby points have nearly the same height), but the walker inches along like a snail. Thousands of steps explore a sliver of the posterior.
• Huge steps — proposals leap into the wasteland and get rejected, so the walker mostly stands still, recording the same value over and over.
• Goldilocks — moderate rejection is the price of real movement. A high acceptance rate is not a good sign by itself.

Computed live in your browser — the same Metropolis math as the notebook (seeded mulberry32; the notebook's centerpiece walker runs in 1-D on the Beta posterior, this page moves it onto a 2-D ridge so you can see exploration). Contours are the target; each colored path is one chain's visit log; the big dots are the walkers now. Changing the step size resets the chains so the diagnostics below stay honest.

Computed live — the trace plot of the same four chains (parameter 1 vs step). At goldilocks the four lines collapse into one overlapping fuzzy band. At 0.02 they become four slow, separate rivers that take ages to meet; at 8 they become flat staircases. Remember these shapes — the clinic below names them.

5 — Trust, but verify: R-hat and ESS on your chains

Everything above worked because we could overlay the exact answer. In real life there is no answer key — the posterior is unknown (that's why we're sampling it!). The standard trick: run several chains from deliberately different starting points. If they all wander to the same distribution and forget where they began, that's strong evidence they found the real posterior; this is mixing — the chain moving freely through the full posterior rather than lingering in one region. If different starts give different answers, the chains are lying to you.

R-hat (\(\hat R\), the Gelman–Rubin statistic) automates the eyeball test: it compares the spread between chains to the spread within chains:

\[ \hat R = \sqrt{\frac{\widehat{\text{var}}^{+}}{W}}, \qquad \widehat{\text{var}}^{+} = \frac{n-1}{n} W + \frac{B}{n}, \]

where \(W\) is the average within-chain variance and \(B\) the between-chain variance. Mixed chains ⇒ between ≈ within ⇒ \(\hat R \approx 1\). The rule of thumb you'll use forever: R-hat below ~1.01 is healthy; above that, do not trust the draws for that parameter.

ESS — effective sample size — is the second auditor. A chain's draws are not independent darts: each position is one small step from the last, so consecutive draws are autocorrelated — correlated with their own recent past. Correlated draws repeat information; ESS is the honest exchange rate: the number of independent draws that would estimate the posterior mean equally well. It's √ESS, not √n, in your error denominator.

Computed live — the autocorrelation of your chain 1, lag by lag: how much does draw t resemble draw t−k? Fast decay = independent-ish draws = high ESS. At step 0.02 the bars barely decay — each draw is nearly a copy of the last. The notebook's pinned 1-D comparison: the tiny-step chain kept an ESS of 19 out of 20,000 draws (0.1%); goldilocks kept 4,205 (21%).

With healthy mixing, the notebook's four dispersed chains (starts 0.05 / 0.35 / 0.65 / 0.95) landed on post-burn-in means of 0.1281 / 0.1262 / 0.1266 / 0.1272 — four corners, one answer — and earned R-hat 1.0009 by hand, 1.0010 by arviz (which uses a stricter rank-normalized split version of the same idea). This is why every PyMC fit runs 2–4 chains by default, and why a single chain is never trustworthy no matter how nice it looks: a stuck chain has no idea it is stuck. Only disagreement between chains exposes the problem.

Make R-hat scream: the two-island trap

R-hat earns its keep when things go wrong. The classic killer is a bimodal posterior — two separate islands of probability (in MMMs this happens when two explanations of the data are both locally plausible). A random-walk sampler can't cross the ocean of near-zero probability between islands: chains that start on different islands stay there, each one locally looking perfectly healthy. Drag the island separation and find the cliff edge:

Computed live in your browser — four seeded chains (two started on each island), 3,000 steps each, classic split R-hat. Small d: one blob, free crossing, R-hat happy. Large d: stranded chains, R-hat explodes. The notebook's measured values: at d = 3 the by-hand split R-hat hit 6.56 (arviz's rank-normalized version: 1.73); at d = 1 arviz gave 1.01; both far past the 1.05 alarm at any real separation.

6 — NUTS: the engine PyMC actually uses (optional deep dive)

Random-walk Metropolis is a blindfolded walker: every proposal is a guess. In 40 dimensions almost every direction is "downhill into the wasteland," so the walker must take minute steps to keep any acceptance at all — autocorrelation explodes and ESS collapses. The curse of dimensionality, back for revenge.

Modern samplers remove the blindfold. The posterior's log-height has a gradient — the local slope, telling you uphill from downhill — and computers get it nearly free (the same automatic differentiation that trains neural networks). Hamiltonian Monte Carlo (HMC) uses it with a physics trick: turn the posterior into a landscape (high probability = low valley), place a puck at the current position, flick it in a random direction, and let it glide — momentum carrying it through long, sweeping, curve-hugging paths instead of dithering. The end of the glide is the proposal: typically far from the start yet still in a high-probability region, so it's usually accepted. Distant, nearly-uncorrelated draws — huge ESS per draw. NUTS (the No-U-Turn Sampler) is HMC with the last knob automated: it lets the puck glide until the path starts to U-turn back toward where it began, then stops. NUTS is the default in PyMC and in mmm_framework (which typically runs it on the fast NumPyro backend — the default, selected on the model config via ModelConfigBuilder().bayesian_numpyro(), not a fit() argument).

Computed live in your browser — the same leapfrog math as the notebook. Same compute budget on the same ridge: 60 random-walk proposals (red, dithering near the start) vs one 60-gradient-step glide (green, sweeping along the valley to a distant, still-probable proposal ★). The notebook's seeded run measured farthest reach 5.2 for the puck vs 1.4 for the walker, with energy drift along the glide |ΔH| = 0.011 — small drift means the distant proposal is almost always accepted.

The notebook then hands the same conversion-rate model to PyMC — same prior, same data, same posterior our 20-line walker sampled — and lets NUTS at it:

Measured — the notebook's seeded comparison table. Read it honestly: a 1-D target is random walk's best case, and NUTS still wins on information per draw — the column that matters. The gap becomes a chasm as dimensions grow: random-walk efficiency collapses with dimension while gradient-guided gliding barely notices. Your 40-parameter MMM in workshop 03 is only feasible because of this.

Divergences: the engine warning light

The table has a column Metropolis can't have: divergences — moments when the gliding puck's simulated physics broke down: the trajectory shot off the landscape and the step had to be abandoned. They happen in regions of extreme curvature, where the glide step is too coarse to follow the surface — a canyon whose walls bend sharper than the puck's stride. Three things for an analyst to burn in:

1. A divergence is a location, not a nuisance. It marks a specific region the sampler tried and failed to enter — your draws may be blind to exactly that region, so the posterior you see can be biased, not just noisy.
2. Don't ignore them; don't just silence them. A handful in tens of thousands of draws may be tolerable; dozens are a model problem. The real fixes change the geometry: gentler/tighter priors that don't carve cliffs, or reparameterization — rewriting the model in equivalent variables whose landscape is smoother (the famous "non-centered" trick for hierarchical models).
3. You mostly inherit the fixes. mmm_framework's default priors (see workshop 01) and its hierarchical structures already encode these reparameterizations — one of several reasons workshop 03 starts from the framework's defaults instead of hand-rolled PyMC.

7 — The trace-plot reading clinic

A trace plot — parameter value vs step number, one line per chain — is the sampler's EKG, and reading one is a 10-second skill that catches most problems before any statistic does. Four patterns cover ~90% of what you'll ever see. All four panels below are real chains, generated in your browser with the notebook's exact recipes (same targets, same step sizes):

Computed live in your browser (seeded) — the notebook's exact recipes. The notebook's measured diagnostics for its versions of these four panels: healthy R-hat 1.00; trending 3.93; stuck 1.74; and the sticky panel passes R-hat but pays in ESS — 74 vs the healthy panel's 1,099. The eyeball and the statistics agree.

8 — The crucial caveat: green lights validate the computation, not the model

Treat today's diagnostics as the pre-flight checklist, not the navigation. The full treatment of this trap — 16 synthetic worlds, 8 silent attribution failures, every one with r-hat ≤ 1.02 and zero divergences — is the pressure-testing scorecard, and its opening guide The Rosy Picture shows two worlds with identical green dashboards and opposite verdicts. And it's why workshop 03 fits your first MMM on a synthetic world whose true answer is known: the only way to learn what a posterior should look like is to practice where truth is checkable.

9 — Glossary