Simulation-based Bayes

Why this matters

The beta-binomial model has a neat closed-form posterior. Most real models do not: a logistic regression with a dozen predictors, a hierarchical model of many hospitals or machines, a pharmacokinetic curve. Their posteriors have no formula, and the normalising constant P(data) is an intractable integral. Markov chain Monte Carlo (MCMC) gets around this by drawing samples from the posterior without ever computing that constant. Tools such as PyMC and Stan, used across data science, biostatistics and actuarial modelling, are built on it.

The key insight

Bayes' rule says

posterior(θ) = likelihood(θ) × prior(θ) / P(data)

P(data) does not depend on θ. So the ratio of the posterior at two points needs only the numerator:

posterior(θ_new) / posterior(θ_old) = [likelihood × prior](θ_new) / [likelihood × prior](θ_old)

If we can compare points, we can build a random walk that spends time in each region in proportion to its posterior probability. The positions it visits are then (dependent) samples from the posterior, and everything from module 2 applies: means, percentiles and probabilities come from counting samples.

The Metropolis algorithm

1. Start at some value θ.
2. Propose a nearby value: θ' = θ + noise, with noise drawn from Normal(0, step).
3. Compute the acceptance ratio r = posterior(θ') / posterior(θ).
4. Accept the move with probability min(1, r): always move uphill; sometimes move downhill.
5. Record the current position (the new one if accepted, otherwise the old one again).
6. Repeat thousands of times.

Work with log posteriors to avoid underflow: accept when random < exp(min(0, log_post(θ') - log_post(θ))). A proposal outside the allowed range (such as p < 0 for a rate) has log posterior -∞, so exp(-∞) = 0 and it is always rejected.

Worked example: a rate with a known answer

A coin (or a pass/fail process) gives 7 successes in 10 trials, with a Beta(2, 2) prior. From the previous lesson we know the exact posterior is Beta(9, 5) with mean 9/14 ≈ 0.643 and SD about 0.124. That lets us check MCMC against the truth.

import math, random

def log_post(p):
    if not 0 < p < 1:
        return -math.inf
    # log prior Beta(2, 2) + log likelihood (7 of 10), constants dropped
    return (1 + 7) * math.log(p) + (1 + 3) * math.log(1 - p)

rng = random.Random(42)
current, current_lp = 0.5, log_post(0.5)
samples, accepted = [], 0
for _ in range(30_000):
    proposal = current + rng.gauss(0, 0.1)
    proposal_lp = log_post(proposal)
    if rng.random() < math.exp(min(0, proposal_lp - current_lp)):
        current, current_lp = proposal, proposal_lp
        accepted += 1
    samples.append(current)

kept = samples[2_000:]                 # discard burn-in
print(sum(kept) / len(kept))           # about 0.64
print(accepted / 30_000)               # acceptance rate, about 0.77

The sample mean lands close to 0.643. Replace log_post with any model's log-likelihood plus log-prior and the same loop works.

Tuning and diagnostics

• Burn-in. Early samples depend on the starting point. Discard the first part of the chain (here 2,000 steps).
• Step size. Too small and nearly every move is accepted but the chain crawls; too large and almost everything is rejected. For one parameter, an acceptance rate of very roughly 40 to 80% is workable; tools tune this automatically.
• Autocorrelation. Consecutive samples are dependent, so 28,000 kept samples carry less information than 28,000 independent ones. Libraries report an effective sample size (ESS).
• Convergence. Run several chains from different starts. If they disagree (the R-hat diagnostic is well above 1.01), the chain has not converged and the results are not trustworthy.
• Trace plots. Plot samples against step number: a healthy chain looks like a fuzzy horizontal band, not a slow drift.

In practice

You will rarely hand-write Metropolis for real work. PyMC (Python) and Stan use Hamiltonian Monte Carlo (NUTS), a smarter proposal that handles many parameters efficiently:

import pymc as pm
with pm.Model():
    p = pm.Beta("p", 2, 2)
    pm.Binomial("y", n=10, p=p, observed=7)
    trace = pm.sample(2000, chains=4, random_seed=1)

Writing the simple version once means you can read their diagnostics and know what can go wrong.

Common mistakes

• Forgetting to repeat the current value on rejection. Recording only accepted moves biases the samples.
• Using raw probabilities instead of logs, which underflow for real datasets.
• No burn-in and one short chain. Always check convergence.
• Treating autocorrelated samples as independent when reporting precision.

Recap

• MCMC samples a posterior using only likelihood × prior, never P(data).
• Metropolis: propose nearby, accept with probability min(1, ratio), record the current value every step.
• Discard burn-in, tune the step size, and check convergence with several chains.
• Production tools (PyMC, Stan) use the same idea with smarter proposals.

Try it

In the Do tab, write the Metropolis acceptance step and check your samples against the exact Beta(9, 5) posterior.