5.2 — Bayesian Regression and Prior/Posterior Predictive Checks

Chapter 5 · Rethinking Statistics with Bayesian Methods

Prof. Xuhu Wan

Section 5.2 · Chapter 5 · Learning Statistics with Python

Bayesian Regression and Prior/Posterior Predictive Checks

Rethinking Statistics with Bayesian Methods

Prof. Xuhu Wan

ISOM, HKUST Business School · 2026 Edition

Bayesian Regression and Prior/Posterior Predictive Checks

A regression is a statement that \(y\) is normal with a mean that moves with \(x\). You will simulate from the prior before touching the data, discover which prior is absurd, then simulate from the posterior to see whether the model can reproduce what it was fitted to.

Regression from the normal-distribution perspective

\[y_i \sim N(\mu_i,\ \sigma^2), \qquad \mu_i = \beta_0 + \beta_1 x_i\]

The notebook’s artificial data: \(\beta_0 = 1\), \(\beta_1 = 2\), \(\sigma = 0.5\).

In the Bayesian regression, which quantities get a prior?

  • \(x\) and \(y\)
  • \(y\) only
  • \(\beta_0\), \(\beta_1\) and \(\sigma\)
  • The residuals \(\varepsilon_i\)

Generate and look

y mean 2.041, sd 0.759. The sd of \(y\) (0.76) exceeds \(\sigma\) (0.5) because part of the spread is the line itself. Keep that distinction: spread of y ≠ noise.

The model in PyMC (Colab only)

In Colab

The notebook fits this with NUTS. We will get the same posterior in closed form on the next slides.

with pm.Model() as linear_regression1:
    std   = pm.HalfCauchy("sigma", beta=10, initval=1)
    beta0 = pm.Normal("intercept", 0, sigma=20)
    beta1 = pm.Normal("slope", 0, sigma=20)
    likelihood = pm.Normal("y", mu=beta0 + beta1 * data["x"], sigma=std,
                           observed=data["y"])
    trace1 = pm.sample(tune=6000, draws=2500, chains=2)

The posterior draws come back as a table: one row per draw, one column per parameter. Every plot in this section is built from that table.

Prior predictive check

Before fitting, draw parameters from the prior and plot the regression lines they imply. If the lines are absurd, the prior is absurd — and you found out for free.

A second data set, wider scale

\(x \in [-10, 10]\), \(y = 0.5 + 3x + N(0, 6^2)\), 200 points.

\(x\) has sd 5.8 and the slope is 3, noise sd 6. Roughly what is the sd of \(y\)? (\(\sqrt{3^2 \cdot 5.8^2 + 6^2}\))

round(np.sqrt(9 * 5.8**2 + 36), 1)

18.4

Prior I: N(0, 10) on both coefficients

Draw 100 lines with \(\beta_0, \beta_1 \sim N(0, 10)\). At \(x = 10\), what range covers 95 % of them?

  • About ±3
  • About ±10
  • About ±30
  • About ±200

Prior II: N(0, 1) and N(0, 1.2)

Prior I: [-173.8, 209.7]. Prior II: [-20.8, 25.2]. The data at \(x = 10\) sit near 30. Prior II covers plausible slopes without covering nonsense — but it puts the true slope 3 at 2.5 prior sd. With 200 points the data will drag it there anyway; with 10 points they might not.

Posterior in closed form

Normal likelihood, normal prior on \(\beta\): the posterior is normal and you can write it down. The formula is the same precision-weighting you saw on the grid, now in matrix form.

Precision adds, means average

\[\Sigma_n = \left(\Sigma_0^{-1} + \tfrac{1}{\sigma^2} X^\top X\right)^{-1}, \qquad \mu_n = \Sigma_n\left(\Sigma_0^{-1}\mu_0 + \tfrac{1}{\sigma^2} X^\top y\right)\]

OLS on (xt, yt) gives slope 3.011. Prior II is \(N(0, 1.2^2)\) on the slope, and the data precision on the slope is about \(\sum x^2/\sigma^2 \approx 6700/32.5 \approx 206\). Will the posterior mean slope be nearer 3.011 or 0?

posterior mean slope ≈ ?

≈ 3.00 — data precision 206 vs prior precision 0.69

Compute it

OLS [0.385, 3.011], posterior mean [0.331, 3.000], posterior sd [0.374, 0.069]. The intercept is shrunk a little toward 0 (its prior sd 1 competes with a data sd 0.4); the slope barely moves. With \(n = 200\) the prior is a rounding error.

Sampling the posterior — one row is one line

The 50 red lines are 50 rows of post_df. Slope credible interval [2.861, 3.139] — the same information a classical CI gives, but readable as a probability.

Posterior predictive check

Now go one step further: for each posterior draw, simulate a whole replicated data set. If the replicas do not look like the data, the model is missing something.

What is one row of the PPC?

For posterior draw \(s\): \(\mu_i^{(s)} = b_0^{(s)} + b_1^{(s)}x_i\), then \(y_i^{\text{rep},(s)} \sim N(\mu_i^{(s)}, \sigma)\) for every \(i\).

yrep has shape (1000, 200). What is yrep[7]?

  • 1000 predictions for observation 7
  • One simulated data set of 200 points, from posterior draw 7
  • The 7th observed value repeated
  • The posterior mean line

Build it and compare

Thirty cyan replicas and the black data share one shape. If the data had a bump the replicas never produce, the normal likelihood would be wrong — this is how the §5.3 outliers will announce themselves.

Turn the picture into a number

A Bayesian p-value: the fraction of replicas whose statistic exceeds the observed statistic. Near 0.5 is good; near 0 or 1 means the model cannot produce what you saw.

p_mean is done. Fix p_max so it compares the maximum of each replica row with the observed maximum yt.max(). Expect about 0.13 — the normal model slightly under-produces the largest value.

What you discovered

  • Regression is a distribution: \(y \mid x \sim N(\beta_0 + \beta_1 x, \sigma^2)\), with priors on \(\beta\) and \(\sigma\).
  • A prior predictive check costs nothing and showed Prior I spraying lines over ±200 where the data live in ±35.
  • With a normal likelihood the posterior is closed-form: precisions add, means are precision-weighted. Slope 3.011 → 3.000.
  • One posterior draw = one regression line; one PPC row = one replicated data set.
  • A Bayesian p-value near 0.5 means the model reproduces that feature; p_max = 0.13 is the first hint that normal tails are thin.

Next: §5.3 — put three outliers in and watch the normal likelihood surrender.

Prof. Xuhu Wan · HKUST ISOM · Learning Statistics with Python