Normal model with exchangeable parameters#

Part 1 · Stage 5 · 🏛️ Hierarchical Models · Lesson 036 of 144 · beginner

◀ Previous · Bayesian analysis of conjugate hierarchical models · Next · Example: parallel experiments in eight schools ▶ · ↑ 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.

The hierarchical normal#

Replace binomial groups with normal ones and the hierarchy becomes fully transparent — every quantity has a closed form, and the mechanics of shrinkage can be read off directly. Each of \(J\) groups supplies an estimate \(\bar{y}_j\) of its own mean \(\theta_j\), with known standard error \(\sigma_j\):

\[\bar{y}_j \mid \theta_j \sim \mathrm{N}(\theta_j, \sigma_j^2), \qquad \theta_j \mid \mu, \tau \sim \mathrm{N}(\mu, \tau^2), \qquad (\mu, \tau) \sim p(\mu, \tau).\]

Here \(\mu\) is the population mean and \(\tau\) the between-group standard deviation: how different the groups really are.

Shrinkage, in closed form#

Conditional on \((\mu, \tau)\), each group’s posterior is the Stage 2 normal update — precisions add, and the mean is precision-weighted:

\[\theta_j \mid \mu, \tau, y \;\sim\; \mathrm{N}(\hat{\theta}_j,\; V_j), \qquad \hat{\theta}_j = \frac{\frac{1}{\sigma_j^2}\bar{y}_j + \frac{1}{\tau^2}\mu} {\frac{1}{\sigma_j^2} + \frac{1}{\tau^2}}, \qquad \frac{1}{V_j} = \frac{1}{\sigma_j^2} + \frac{1}{\tau^2}.\]

Define the shrinkage factor \(B_j = \sigma_j^2 / (\sigma_j^2 + \tau^2)\); then \(\hat{\theta}_j = (1 - B_j)\, \bar{y}_j + B_j\, \mu\). A noisy group (large \(\sigma_j\)) has \(B_j\) near 1 and is pulled almost entirely to the population mean; a precise group keeps its own estimate. And crucially, \(\tau\) is not chosen — it is inferred, so the data decide how much pooling is warranted.

The whole model in code#

import pymc as pm
with pm.Model():
    mu  = pm.Normal("mu", 0, 10)
    tau = pm.HalfNormal("tau", 10)               # weakly informative; never inverse-gamma(eps,eps)
    theta = pm.Normal("theta", mu, tau, shape=J)
    pm.Normal("y", theta, sigma=sigma_j, observed=ybar)   # sigma_j known
    idata = pm.sample(target_accept=0.95)

Two warnings#

When \(\tau\) is near zero, the posterior sits at a boundary — the counterexample from Stage 4 — and the geometry becomes a funnel that samplers negotiate badly; the standard repair is the non-centred parameterisation, \(\theta_j = \mu + \tau \eta_j\) with \(\eta_j \sim \mathrm{N}(0,1)\). And the prior on \(\tau\) matters, especially with few groups: with \(J = 8\), a careless inverse-gamma can dominate. Both issues are met head-on in the eight-schools example that follows.

Hint