11  Lab 10: Priors, Posteriors, and Checking the Sampler

Run this one live: the executable lab.

11.1 What we are doing

Every model so far returned a number. Today’s returns a distribution, and the work shifts accordingly: less about fitting, more about deciding what you believe beforehand and checking whether the machinery actually worked.

The question is small on purpose so the workflow stays visible: what fraction of weeks are high-SNAP weeks at CA_1 FOODS? You will state a prior, check it before looking at the data, fit the model, verify the sampler, and then ask the posterior questions no point estimate can answer.

Suggested time: ~35 minutes · You will need: pandas, numpy, matplotlib, pymc, arviz, and data/processed/m5_weekly.csv. Each fit takes about ten seconds.

---

11.2 Setup

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import pymc as pm
import arviz as az
import warnings, logging
warnings.filterwarnings('ignore')
logging.getLogger('pymc').setLevel(logging.ERROR)

from pathlib import Path

def data_path(name):
    for base in (Path.cwd(), *Path.cwd().parents):
        p = base / "data" / "processed" / name
        if p.exists():
            return p
    raise FileNotFoundError(
        f"Could not find data/processed/{name}. "
        "Run `python scripts/prep_m5.py` from the repository root first.")

weekly = pd.read_csv(data_path("m5_weekly.csv"), parse_dates=["week_start"])
ca1 = weekly[(weekly.store_id == "CA_1") & (weekly.cat_id == "FOODS")]

N_WEEKS = len(ca1)
N_HIGH  = int((ca1.snap_days >= 5).sum())
print(f"CA_1 FOODS: {N_HIGH} high-SNAP weeks out of {N_WEEKS} "
      f"({100*N_HIGH/N_WEEKS:.0f}% observed)")

ArviZ changed its default interval

az.summary(idata) now reports an 89% equal-tailed interval (eti89_lb, eti89_ub). This course uses the highest density interval, which is the shortest interval containing a given probability mass.

Always ask for it: az.summary(idata, ci_kind="hdi", ci_prob=0.94).

For a symmetric posterior the two nearly coincide. For a skewed one they do not, and the HDI is the honest summary.

---

Walkthrough

CA_1 FOODS: 72 high-SNAP weeks out of 277 (26% observed)

Seventy-two of two hundred seventy-seven. Twenty-six percent. Hold that number — the whole lab is about how different beliefs collide with it.

Read this callout, because it will bite you in Homework 6 otherwise. ArviZ changed its default interval. az.summary now reports an 89% equal-tailed interval. This course uses the highest density interval — the shortest interval containing the mass — so you must ask for it every time: ci_kind="hdi", ci_prob=0.94.

For a symmetric posterior the two nearly coincide; here the equal-tailed 94% interval is [0.2155, 0.3141] against an HDI of [0.2144, 0.3129], so it hardly matters. For a skewed posterior they diverge, and the HDI is the honest summary. Get in the habit now, while the stakes are low.

11.3 Step 1 — Check the prior before you use it

We will model the count of high-SNAP weeks as Beta-Binomial. Start with Beta(2, 2) — weakly informative, centered at 0.5.

Before fitting anything, ask what that prior alone implies about the data.

with pm.Model() as model:
    theta = pm.Beta("theta", alpha=2, beta=2)
    obs   = pm.Binomial("obs", n=N_WEEKS, p=theta, observed=N_HIGH)
    prior_pred = pm.sample_prior_predictive(500, random_seed=42)

sim = prior_pred.prior_predictive["obs"].values.ravel()
fig, ax = plt.subplots(figsize=(8, 3))
ax.hist(sim, bins=40)
ax.axvline(N_HIGH, color="crimson", lw=2, label=f"observed = {N_HIGH}")
ax.set_xlabel("simulated high-SNAP weeks out of 277"); ax.legend()
plt.tight_layout()
print(f"prior predictive spans {sim.min()} to {sim.max()} weeks")

Read the plot. The prior thinks almost any count is plausible — from nearly zero to nearly every week. The observed value sits comfortably inside, so the prior is not wrong.

But is it good? SNAP is disbursed on a published calendar. We are not genuinely ignorant here, and a prior that entertains “every week is a high-SNAP week” is entertaining something we know is false.

---

Walkthrough

Here is the step almost everyone skips, and it costs about four seconds.

Before fitting anything, we ask what the prior alone implies about the data. Beta(2,2) — weakly informative, centered at one half. Draw five hundred values of θ from it, simulate a count of high-SNAP weeks from each, and look at what comes out.

prior predictive spans 8 to 276 weeks

Eight to two hundred seventy-six. Out of 277. The prior considers essentially every possible answer plausible — from almost never to almost always.

The observed value, 72, sits comfortably inside that. So the prior is not wrong. Nothing here is going to blow up.

But look at what it is willing to entertain. “Two hundred seventy-six of two hundred seventy-seven weeks are high-SNAP weeks” is not an uncertain belief, it is a false one. SNAP is disbursed on a published calendar. We are not genuinely ignorant about this rate, and Beta(2,2) pretends we are.

That is the difference between a prior that is defensible and one that is lazy, and the prior predictive check is what makes the difference visible. It costs one function call, and it is the single highest-value habit in this lecture.

11.4 Step 2 — Fit it

with model:
    idata = pm.sample(2000, tune=1000, chains=4, progressbar=False, random_seed=42)

summary = az.summary(idata, ci_kind="hdi", ci_prob=0.94)
print(summary.round(4).to_string())

The posterior mean sits between the prior mean of 0.5 and the observed proportion — much closer to the data, because 277 observations carry far more weight than a Beta(2, 2) prior.

---

Walkthrough

         mean     sd hdi94_lb hdi94_ub ess_bulk ess_tail r_hat
theta  0.2628  0.026     0.22     0.31     3686     5175  1.00

Posterior mean 0.2628.

Now put three numbers side by side, because their arrangement is the lesson. The prior mean was 0.5. The observed proportion is 0.2599. The posterior mean is 0.2628.

The posterior sits between them, as it must — but look at where between them. It is three thousandths from the data and nearly a quarter away from the prior. Two hundred seventy-seven observations against a Beta(2,2) prior is not a fair fight, and the data wins almost totally.

That is what “the likelihood dominates” means, concretely, and it is the answer to half the objections people raise about Bayesian methods before they have run one.

The 94% HDI reads [0.22, 0.31] here — two decimals, because that is az.summary’s default rounding. Step 4 prints the same interval more precisely. Same numbers, different display.

11.5 Step 3 — Did the sampler work?

MCMC always returns samples. Whether they came from the posterior is a separate question, and one you are expected to answer every time.

n_div = int(idata.sample_stats.diverging.sum())
print(f"R-hat        {summary.r_hat.max():.3f}    (want < 1.01)")
print(f"ESS (bulk)   {summary.ess_bulk.min():.0f}      (want > 400)")
print(f"divergences  {n_div}        (want 0)")

az.plot_trace(idata, var_names=["theta"]);
plt.tight_layout()

Read the trace plot. On the right, four chains overlapping in a band with no drift and no stuck stretches — that is what convergence looks like. On the left, four densities lying on top of one another.

If the chains had separated, \(\hat{R}\) would exceed 1.01 and nothing downstream would mean anything.

---

Walkthrough

R-hat        1.001    (want < 1.01)
ESS (bulk)   3687      (want > 400)
divergences  0        (want 0)

MCMC always returns samples. It never fails with an error because the answer is wrong. Whether those samples came from the posterior you asked for is a separate question, and you are expected to answer it every single time before reading any result.

By now that sentence should sound familiar. It is the same shape as zero_grad(), as the transposed axis, as the missing positional encoding — the code runs, the output is plausible, and nothing tells you. Bayesian workflow just has better instruments for catching it.

Three numbers. R-hat at 1.001, comfortably under 1.01 — the four chains agree with each other. Effective sample size at 3,687, far above the 400 rule of thumb — the samples are not so autocorrelated that you have less information than it looks. Divergences at zero — the sampler never had to give up on a step.

And the picture. On the right, four chains overlapping in a flat band, no drift, no stuck stretches — that is what convergence looks like, and it is worth memorizing the shape so a bad one jumps out. On the left, four densities lying on top of each other.

If those chains had separated, R-hat would exceed 1.01 and nothing downstream would mean anything — not the mean, not the interval, not the probabilities in the next step.

11.6 Step 4 — Ask the posterior a question

The posterior is a set of samples, so any question about \(\theta\) is answered by counting.

theta = idata.posterior["theta"].values.ravel()

print(f"P(theta > 0.30)          = {(theta > 0.30).mean():.3f}")
print(f"P(0.24 < theta < 0.28)   = {((theta > 0.24) & (theta < 0.28)).mean():.3f}")
print(f"94% HDI                  = {np.round(az.hdi(theta, prob=0.94), 4)}")

Say the first line out loud as a sentence a manager would understand. Then notice you cannot say anything of that form from a confidence interval — a frequentist 94% interval is a statement about the procedure, not about \(\theta\).

---

Walkthrough

The posterior is a set of samples. So any question about θ is answered by counting, which is the part that feels almost too easy.

P(theta > 0.30)          = 0.080
P(0.24 < theta < 0.28)   = 0.558
94% HDI                  = [0.2132 0.3098]

Try it

— (replaces “say it out loud as a sentence a manager would understand”)

Pause and write the first line as a sentence you could say in a meeting. Then try to write the same sentence from a confidence interval.

The first one is easy: “There is an eight percent chance that more than thirty percent of weeks are high-SNAP weeks.” A claim about the world, in the units of the question.

The second is impossible, and this is the part worth being precise about rather than glib. A frequentist 94% interval is a statement about the procedure — that intervals built this way cover the true value 94% of the time. It is not a statement about θ, and “there is a 94% chance θ is in this interval” is exactly the misreading every textbook warns against and every practitioner makes anyway.

The Bayesian version means what people already thought the frequentist version meant. That is the practical appeal, and it is bought with the prior.

11.7 Step 5 — How much did the prior matter?

The honest objection to Bayesian work is that the prior is arbitrary. Test how much it moves the answer.

rows = []
for a, b, label in [(2, 2, "Beta(2,2) — weak, centered 0.5"),
                    (5, 15, "Beta(5,15) — informed, centered 0.25"),
                    (20, 2, "Beta(20,2) — badly wrong, centered 0.91")]:
    with pm.Model():
        th = pm.Beta("theta", alpha=a, beta=b)
        pm.Binomial("obs", n=N_WEEKS, p=th, observed=N_HIGH)
        i2 = pm.sample(1000, tune=1000, chains=2, progressbar=False, random_seed=42)
    s = az.summary(i2, ci_kind="hdi", ci_prob=0.94)
    rows.append({"prior": label, "prior_mean": round(a/(a+b), 3),
                 "posterior_mean": round(float(s["mean"].iloc[0]), 3)})
print(pd.DataFrame(rows).to_string(index=False))

With enough data the prior stops mattering

Three priors, including one centered on 0.91 — badly, deliberately wrong. All three posteriors land close together, because 277 observations dominate.

That is the answer to the arbitrariness objection, and it comes with a condition attached: it holds when you have enough data. With 8 observations rather than 277, the third prior would still be visible in the answer.

Which is exactly the situation Lecture 11 Part 2 is about.

Your turn: re-run the badly-wrong prior against only the first 10 weeks of data and report how far the posterior moves.

# ---- your turn: uncomment and fill in ----
# n_small, k_small = 10, int((ca1.snap_days >= 5).head(10).sum())
# with pm.Model():
#     th = pm.Beta("theta", alpha=___, beta=___)      # <- the badly-wrong prior
#     pm.Binomial("obs", n=n_small, p=th, observed=k_small)
#     i3 = pm.sample(1000, tune=1000, chains=2, progressbar=False, random_seed=42)
# print(az.summary(i3, ci_kind="hdi", ci_prob=0.94).round(3).to_string())

---

Walkthrough

Now the honest objection, the one you should expect from any audience: the prior is arbitrary, so the answer is arbitrary. Test it rather than argue about it.

Three priors. Weak and centered at 0.5. Informed and centered at 0.25 — roughly right. And one centered at 0.91, which is not a mild disagreement, it is a belief that nine weeks in ten are high-SNAP weeks. Deliberately, badly wrong.

                                 prior  prior_mean  posterior_mean
         Beta(2,2) — weak, centered 0.5       0.500           0.265
   Beta(5,15) — informed, centered 0.25       0.250           0.259
Beta(20,2) — badly wrong, centered 0.91       0.909           0.309

The prior means span 0.25 to 0.91 — a range of 0.66. The posterior means span 0.259 to 0.309 — a range of 0.05. The data compressed the disagreement by a factor of thirteen.

Now let me be more careful than the word “close” allows, because this is the step where it is tempting to oversell.

The badly-wrong prior did not land in the same place as the others. It landed at 0.309, and the weak prior’s 94% interval runs up to 0.313. So the wrong prior’s answer sits just inside the upper edge of what the honest analysis considered plausible — inside the range, but at the very boundary of it. And its own interval, [0.258, 0.358], barely contains the observed proportion of 0.260 at its lower edge.

So the correct claim is not “the prior did not matter.” It is: with 277 observations, even an absurd prior gets dragged to the edge of the right answer rather than away from it. That is a strong result and it does not need embellishment.

And it comes with a condition, which is the sentence to land hardest on. It holds when you have enough data. With eight observations instead of 277, that third prior would still be sitting in your answer, plainly visible.

Try it

That is your turn. Re-run the badly-wrong prior against only the first ten weeks and report how far the posterior moves. Fill in the two blanks — they are the alpha and beta of the Beta(20,2).

And that situation — a handful of observations per unit, and a prior that therefore still matters — is not a contrived exercise. It is Lecture 11 Part 2, where ten series have very different amounts of data and the thin ones borrow strength from the thick ones. Everything you just did by hand, the hierarchy will do automatically.

11.8 Before you leave

Discuss (the prior): Your Step 5 table says a deliberately wrong prior barely changed the answer. Under what circumstances would you not be reassured by that, and what would you check?

Discuss (the check): You verified \(\hat{R}\), ESS and divergences before reading any result. Name the equivalent check you have been doing all semester for a non-Bayesian model, and say which is easier to skip.

Discuss (for Homework 6): Part 1 asks you to run a prior predictive check and propose a better prior than Beta(2, 2) for a rate you know is roughly a quarter. What alpha and beta would you choose, and how would you defend the strength of it rather than just the center?

---

Solutions posted to Canvas after class.

Walkthrough

Try it

Three questions for the board.

The prior. Your Step 5 table says a deliberately wrong prior barely changed the answer. Under what circumstances would you not be reassured by that? Name two, and be specific — one about sample size, one about what the prior is a prior over.

A prior on a variance, or on a hierarchical scale parameter, can dominate a posterior with far more data than 277 observations. Lecture 11 Part 2 has an example.*

The check. You verified R-hat, ESS and divergences before reading any result. Name the equivalent check you have been running all semester on non-Bayesian models — and say which of the two is easier to skip.

is easier to skip, because MCMC hands you a tidy summary table that looks authoritative whether or not it converged.*

For Homework 6. Part 1 asks you to run a prior predictive check and propose something better than Beta(2,2) for a rate you know is around a quarter. Choose your alpha and beta — and defend the strength, not just the center. Beta(5,15) and Beta(50,150) have the same mean and say very different things about how sure you are.

Solutions for the coded parts go up on Canvas after the deadline. Next week: the same machinery pointed at a time series.

11.9 What your run should produce

Your numbers will differ where a model is stochastic; these are the values this lab produces on the pinned environment.

random_seed=42 throughout; ~10 s per fit, three more in Step 5. Conjugate, so every value is verifiable in closed form — the exact column below is Beta(2+72, 2+205) etc.

Quantity	MCMC	Exact
Data	72 high-SNAP weeks / 277 (0.2599)	—
Prior predictive span	8 → 276 weeks	—
Posterior mean	0.2628	0.2633
Posterior sd	0.026	—
94% HDI	[0.2132, 0.3098]	[0.2144, 0.3129]
94% equal-tailed (for contrast)	—	[0.2155, 0.3141]
R-hat / ESS bulk / divergences	1.001 · 3,687 · 0	—
P(θ > 0.30)	0.080	0.084
P(0.24 < θ < 0.28)	0.558	0.553
Beta(2,2) → posterior mean	0.265	0.2633 · HDI [0.2144, 0.3129]
Beta(5,15) → posterior mean	0.259	0.2593 · HDI [0.2119, 0.3072]
Beta(20,2) → posterior mean	0.309	0.3077 · HDI [0.2579, 0.3580]
Prior-mean spread → posterior spread	0.66 → 0.05 (13× compression)	—