Example: hierarchical normal model#

Part 3 · Stage 9 · ⛓️ MCMC: Gibbs, Metropolis & HMC · Lesson 074 of 144 · intermediate

◀ Previous · Eﬀective number of simulation draws · Next · Eﬃcient Gibbs samplers ▶ · ↑ 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.

Gibbs, worked through#

The hierarchical normal model of Stage 5 has closed-form conditionals at every level, which makes it the canonical Gibbs example — and a useful place to watch a sampler misbehave. Take \(J\) groups, \(n_j\) observations each:

\[y_{ij} \sim \mathrm{N}(\theta_j, \sigma^2), \qquad \theta_j \sim \mathrm{N}(\mu, \tau^2), \qquad p(\mu, \log \sigma, \log \tau) \propto 1 .\]

The unknowns are the group means \(\theta\), the population mean \(\mu\), and the two variances.

The four conditionals#

Each is a Stage 2 or Stage 3 result, reused:

\[\theta_j \mid \cdot \sim \mathrm{N}\!\left( \frac{\frac{n_j}{\sigma^2}\bar{y}_{\cdot j} + \frac{1}{\tau^2}\mu} {\frac{n_j}{\sigma^2} + \frac{1}{\tau^2}}, \; \left(\frac{n_j}{\sigma^2} + \frac{1}{\tau^2}\right)^{-1}\right), \qquad \mu \mid \cdot \sim \mathrm{N}\!\left(\bar{\theta}, \; \frac{\tau^2}{J}\right),\]

with \(\sigma^2\) and \(\tau^2\) drawn from scaled inverse-\(\chi^2\) distributions built from the within- and between-group sums of squares. Cycle through them and the chain converges to the joint posterior — no tuning, no rejections.

import numpy as np
from scipy import stats
for t in range(n_iter):
    V = 1.0 / (n_j / sigma**2 + 1.0 / tau**2)                  # theta | mu, sigma, tau, y
    m = V * (n_j * ybar_j / sigma**2 + mu / tau**2)
    theta = stats.norm(m, np.sqrt(V)).rvs()
    mu = stats.norm(theta.mean(), tau / np.sqrt(J)).rvs()      # mu | theta, tau
    ss_w = ((y - theta[group]) ** 2).sum()                     # sigma^2 | theta, y
    sigma = np.sqrt(ss_w / stats.chi2(n).rvs())
    ss_b = ((theta - mu) ** 2).sum()                           # tau^2 | theta, mu
    tau = np.sqrt(ss_b / stats.chi2(J - 1).rvs())

Watch it stick#

Now look at what the chain does when \(\tau\) wanders near zero. The conditional for \(\theta_j\) collapses onto \(\mu\); the conditional for \(\tau^2\) is then built from \(\theta_j\) that are all nearly equal, so it stays small. The two conditionals trap each other, and the sampler crawls through the neck of the funnel from eight schools. Gibbs does not diverge or warn — it simply mixes so slowly that the low-\(\tau\) region is under-visited, and \(\hat{R}\) may look fine on a short run.

The lesson#

Closed-form conditionals guarantee correctness in the limit, not efficiency in practice. The diagnosis is posterior correlation between levels, and the two cures are the subject of the next lesson: reparameterise so the levels decouple, or update them jointly.

Hint