Example: hierarchical normal model (continued)#

Part 3 · Stage 10 · 🎛️ Modal & Variational Approximation · Lesson 086 of 144 · intermediate

◀ Previous · Conditional and marginal posterior approximations · Next · Variational inference ▶ · ↑ 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 same model, approximated#

Return to the hierarchical normal model that Stage 9 sampled with Gibbs, and apply the machinery of this stage instead. It is a fair test: the answers are known, so every approximation can be scored.

\[y_{ij} \sim \mathrm{N}(\theta_j, \sigma^2), \qquad \theta_j \sim \mathrm{N}(\mu, \tau^2), \qquad \phi = (\mu, \sigma, \tau).\]

Why the joint mode fails#

Maximise \(p(\theta, \phi \mid y)\) jointly and the optimiser walks straight to the degenerate point: every \(\theta_j\) equal to \(\mu\), and \(\tau = 0\). The joint density there is unbounded — the likelihood of the group parameters concentrates without limit as the population variance shrinks. The joint mode is not merely a poor summary; it does not exist as a finite maximum.

Marginal mode, via EM#

Integrate the \(\theta_j\) out and maximise the marginal \(p(\phi \mid y)\) instead. EM does this without doing the integral: the E-step supplies conditional means and variances of the \(\theta_j\), and the M-step’s update for \(\tau^2\) includes those variances, so it cannot collapse to zero. The degeneracy disappears the moment the nuisance parameters are averaged rather than maximised.

import numpy as np
for _ in range(200):
    V = 1 / (n_j / sigma**2 + 1 / tau**2)                  # E-step
    Etheta = V * (n_j * ybar_j / sigma**2 + mu / tau**2)
    mu = Etheta.mean()                                     # M-step
    tau = np.sqrt(np.mean((Etheta - mu) ** 2 + V))         # + V prevents collapse
# then: Laplace around (mu, log sigma, log tau) for approximate uncertainty

Then the conditional#

Having a marginal for \(\phi\) — from EM plus curvature, or from a grid — the previous lesson’s recipe completes the picture: draw \(\phi\), then draw each \(\theta_j\) exactly from its normal conditional. The group-level answers inherit the hyperparameter uncertainty, which is what distinguishes this from empirical Bayes.

What the comparison teaches#

Against long-run MCMC, the modal approximation is accurate for \(\mu\) and for the \(\theta_j\) in data-rich groups, and least accurate for \(\tau\) — the parameter whose posterior is skewed, bounded below, and often heaped near zero. That is the general pattern: approximations fail on the variance parameters of hierarchies, which are precisely the parameters that decide how much pooling occurs. Use modal methods for speed and starting values; check the conclusions that hinge on \(\tau\) against a sampler.

Hint