Posterior predictive checking#

Part 2 · Stage 6 · 🔍 Model Checking & Comparison · Lesson 042 of 144 · intermediate

◀ Previous · Do the Inferences from the Model Make Sense? · Next · Graphical posterior predictive checks ▶ · ↑ Section

Important

✨ AI-generated content. This page was written with the assistance of an AI language model and is provided as a learning aid. Despite careful review, it may still contain mistakes, omissions, or out-of-date information. Whether you are new to the topic, a team lead, or a senior practitioner, treat it as a starting point rather than an authoritative reference: read it critically and independently verify anything you act on (code, commands, figures, and factual claims) against official documentation and primary sources before relying on it.

Let the model generate data#

The central technique of Part II is disarmingly simple: if the model fits, data simulated from it should look like the data you actually saw. Simulate many replicated datasets from the fitted model, compare them to the real one, and any systematic difference is a failure of the model.

The posterior predictive distribution#

Replicated data \(y^{\text{rep}}\) are drawn from the posterior predictive distribution — the likelihood averaged over posterior uncertainty in the parameters:

\[p(y^{\text{rep}} \mid y) = \int p(y^{\text{rep}} \mid \theta) \; p(\theta \mid y) \; d\theta .\]

In simulation this is one line per draw: for each posterior draw \(\theta^{(s)}\), generate a full dataset \(y^{\text{rep}(s)} \sim p(y \mid \theta^{(s)})\). Note what is averaged in — parameter uncertainty and sampling variability, exactly the two sources present in the real data.

Test quantities#

You cannot compare whole datasets by eye at scale, so choose a test quantity \(T(y, \theta)\) summarising the feature you care about: the maximum, the number of zeros, a lag-1 autocorrelation, the variance-to-mean ratio. Unlike a classical test statistic, \(T\) may depend on the parameters. The posterior predictive \(p\)-value is the probability that the replicated data are more extreme than what you observed:

\[p_B = \Pr\bigl(T(y^{\text{rep}}, \theta) \ge T(y, \theta) \;\bigm|\; y\bigr) \;\approx\; \frac{1}{S} \sum_{s=1}^{S} \mathbf{1}\bigl\{T(y^{\text{rep}(s)}, \theta^{(s)}) \ge T(y, \theta^{(s)})\bigr\} .\]

Values near 0 or 1 mean the observed data are atypical of the model in that respect.

import numpy as np, pymc as pm
with model:
    ppc = pm.sample_posterior_predictive(idata)
yrep = ppc.posterior_predictive["y"].values.reshape(-1, len(y))

T = lambda d: d.max()                                # any feature you care about
p_B = (np.array([T(r) for r in yrep]) >= T(y)).mean()
p_B                                                  # ≈ 0.5: unremarkable; ≈ 0 or 1: misfit

Read it as a diagnostic, not a verdict#

Two honest caveats. The check uses the data twice — to fit and to test — so \(p_B\) is conservative, tending to concentrate near 0.5; it is not calibrated like a frequentist \(p\)-value and is not “the probability the model is true”. And a check can only find what its test quantity looks for: choosing \(T\) = the sample mean will almost always pass, because the model was fitted to match the mean. Choose test quantities the model does not automatically reproduce, and ones that matter for your purpose. The point is never to accept or reject, but to learn where the model misses.

Hint