1. Model setup: hierarchical normal model
We consider the hierarchical normal model from earlier.
For groups $j = 1, \dots, J$ and observations $i = 1, \dots, n_j$:
- Data model (within-group observations):
- $y_{ij} \mid \theta_j, \sigma^2 \sim \text{Normal}(\theta_j, \sigma^2)$
- Group-level model (between-group variation):
- $\theta_j \mid \mu, \tau^2 \sim \text{Normal}(\mu, \tau^2)$
Parameters:
- $\theta = (\theta_1, \dots, \theta_J)$: group-specific means
- $\mu$: overall mean across groups
- $\sigma$: within-group standard deviation
- $\tau$: between-group standard deviation
For numerical stability and unconstrained optimization, we work with $\log \sigma$ and $\log \tau$, rather than σ and τ directly.
Goal:
Obtain approximate Bayesian inference using mode-based approximations, by
- Finding the joint posterior mode of $(\theta, \mu, \log \sigma, \log \tau)$ and
- Then focusing on the marginal posterior mode of $(\mu, \log \sigma, \log \tau)$, with a normal approximation around that mode.
Given $(\mu, \log \sigma, \log \tau)$, the conditional posterior for each $\theta_j$ is normal and independent across j.
2. Crude initial estimates
We first compute simple, non-Bayesian estimates to use as starting values.
For each group $j$:
- Let $\bar{y}_{\cdot j}$ (often written $y_{\cdot j}$) be the mean of group j:
- $\bar{y}_{\cdot j} = \frac{1}{n_j}\sum_{i=1}^{n_j} y_{ij}$
- Let $s_j^2$ be the within-group sample variance:
- $s_j^2 = \frac{1}{n_j – 1} \sum_{i=1}^{n_j} (y_{ij} – \bar{y}_{\cdot j})^2$
Then define crude estimates:
- Group means: $\hat{\theta}_j^{(0)} = \bar{y}_{\cdot j}$
- Within-group variance σ²: average the within-group variances:
- $\hat{\sigma}^2{}^{(0)} = \frac{1}{J} \sum_{j=1}^J s_j^2$
- Overall mean μ: mean of the $\hat{\theta}_j^{(0)}$:
- $\hat{\mu}^{(0)} = \frac{1}{J} \sum_{j=1}^J \hat{\theta}_j^{(0)}$
- Between-group variance τ²: variance of the $\hat{\theta}_j^{(0)}$:
- $\hat{\tau}^2{}^{(0)} = \text{sample variance of } \{\hat{\theta}_1^{(0)}, \dots, \hat{\theta}_J^{(0)}\}$
In the coagulation example used in the book, these crude estimates appear in the first column of Table. They are already quite close to the final joint mode because there is very little shrinkage (within-group σ² is large relative to τ²).
3. Conditional maximization: finding the joint posterior mode
We now treat the joint posterior $p(\theta, \mu, \log \sigma, \log \tau \mid y)$
and use conditional maximization (stepwise ascent) to find a local joint mode.
Because this hierarchical normal model is conjugate, all conditional posterior distributions have closed forms:
- $p(\theta_j \mid \mu, \sigma, \tau, y)$: Normal
- $p(\mu \mid \theta, \sigma, \tau, y)$: Normal
- $p(\sigma^2 \mid \theta, \mu, \tau, y)$: scaled inverse–χ²
- $p(\tau^2 \mid \theta, \mu, \sigma, y)$: scaled inverse–χ²
We update each parameter by its conditional posterior mode, one parameter (or block) at a time.
3.1 Step 1: conditional mode for each $θ_j$
From the Gibbs-sampling derivations, the conditional posterior of $θ_j$ is:
$\theta_j \mid \mu, \sigma^2, \tau^2, y_j \sim \text{Normal}\left(\hat{\theta}_j, V_{\theta_j}\right)$
where:
- Conditional posterior mean:
- $\hat{\theta}_j = V_{\theta_j} \left(\frac{n_j \bar{y}_{\cdot j}}{\sigma^2} + \frac{\mu}{\tau^2}\right)$
- Conditional posterior variance:
- $V_{\theta_j} = \left(\frac{n_j}{\sigma^2} + \frac{1}{\tau^2}\right)^{-1}$
Since the conditional posterior is normal, the mode equals the mean. So to maximize the joint posterior with respect to $θ_j$, we simply replace $θ_j$ by $\hat{\theta}_j$.
Do this for all j: $j = 1, \dots, J$
3.2 Step 2: conditional mode for μ
Given θ, σ, τ and y, the conditional posterior for μ is also normal:
$\mu \mid \theta, \sigma, \tau, y \sim \text{Normal}\left(\hat{\mu}, V_\mu\right)$
with:
- Posterior mean (and mode): $\hat{\mu} = \frac{1}{J} \sum_{j=1}^J \theta_j$ (assuming a flat or weak prior on μ.
So in conditional maximization, we update: $\mu \leftarrow \hat{\mu}.$
3.3 Step 3: conditional mode for log σ
The conditional posterior of σ² (given θ, μ, τ, y) is scaled inverse–χ².
If we maximize the posterior with respect to σ², we get:
- Mode for σ²:
- $\hat{\sigma}^2_{\text{mode}} = \frac{1}{n} \sum_{j=1}^J \sum_{i=1}^{n_j} (y_{ij} – \theta_j)^2,$
- up to the usual small factors depending on the prior and degrees of freedom.
However, we are working with log σ, not σ². Because of the Jacobian when transforming from σ² to log σ, the mode of log σ is not exactly log of the σ²-mode multiplied by some n/(n+2) factor; instead, the transformation changes the location of the mode slightly. The result (as stated in the text) is:
- Replace current log σ by $\log \hat{\sigma}$ where $\hat{\sigma}$ is defined, i.e. the mode appropriate for log σ given the prior and the Jacobian.
The key idea:
We update log σ to the conditional mode of its own conditional posterior $p(\log \sigma \mid \theta, \mu, \tau, y)$.
3.4 Step 4: conditional mode for log τ
Similarly, the conditional posterior for τ² (given θ, μ, σ, y) is scaled inverse–χ². We compute:
- Mode for τ² (similar structure to σ² case),
- Then adjust for the transformation to log τ, giving a mode $\hat{\tau}$ for τ.
We then update: $\log \tau \leftarrow \log \hat{\tau}$
with $\hat{\tau}$ defined.
3.5 Convergence of conditional maximization
We repeat 1–4 until the parameters stop changing (or change only very slightly). Each step increases the joint posterior density, so the sequence converges to a local joint mode of $p(\theta, \mu, \log \sigma, \log \tau \mid y).$
In the coagulation example:
- Only about 3 iterations are needed to reach approximate convergence.
- The resulting joint mode is very close to the crude initial estimates (because the shrinkage factors $\sigma^2 / (n_j \tau^2 + \sigma^2)$ are near 1: little pooling).
There is also another mode at the boundary of parameter space (e.g., τ → 0), but it has very little posterior mass, so it is not practically important.
4. Why the joint mode may not be a good summary
When J is large but each $n_j$ is small, the joint posterior over $(\theta, \mu, \sigma, \tau)$ can:
- Be highly skewed,
- Have strong dependencies,
- Have modes that don’t represent the typical behavior of the whole posterior mass.
The joint mode is a single point in a high-dimensional space and often does not represent what Bayesian inference cares about most (posterior averages, intervals, tail behavior, etc.).
A better approach:
Focus on a lower-dimensional subset of parameters and approximate their marginal posterior.
So we partition:
- $\gamma = (\theta_1, \dots, \theta_J) = \theta$ (group-level parameters),
- $\phi = (\mu, \log \sigma, \log \tau)$ (hyperparameters).
We then write:
$p(\theta, \mu, \log \sigma, \log \tau \mid y) = p(\theta \mid \mu, \log \sigma, \log \tau, y)\, p(\mu, \log \sigma, \log \tau \mid y)$
Our main target becomes the 3-dimensional marginal:
$p(\mu, \log \sigma, \log \tau \mid y)$
which is easier to approximate with normal (or t) distributions.
5. Exact conditional p(θ | μ, σ, τ, y) and the marginal p(μ, log σ, log τ | y)
From the Gibbs sampler derivation, the conditional posterior of θ given μ, σ, τ, y is:
$p(\theta \mid \mu, \sigma, \tau, y) = \prod_{j=1}^J \text{Normal}\big(\theta_j \mid \hat{\theta}_j, V_{\theta_j}\big)$
with $\hat{\theta}_j, V_{\theta_j}$ as before (functions of μ, σ, τ, and y).
So this conditional is exact, not just approximate.
5.1 Using the ratio identity
We can use the identity from earlier:
$p(\mu, \log \sigma, \log \tau \mid y) = \frac{p(\theta, \mu, \log \sigma, \log \tau \mid y)} {p(\theta \mid \mu, \log \sigma, \log \tau, y)}$
Because the denominator is the exact conditional, this identity holds for any θ. To simplify, we choose θ to be the conditional means:
$\theta_j = \hat{\theta}_j(\mu, \sigma, \tau) \quad\text{for all } j$
Plugging in:
- The joint posterior (up to a constant) is:
- $p(\theta, \mu, \log \sigma, \log \tau \mid y) \propto \tau \prod_{j=1}^J \text{Normal}(\theta_j \mid \mu, \tau^2) \prod_{j=1}^J \prod_{i=1}^{n_j} \text{Normal}(y_{ij} \mid \theta_j, \sigma^2)$
- (ignoring constants not depending on parameters).
- The conditional posterior:
- $p(\theta \mid \mu, \log \sigma, \log \tau, y) = \prod_{j=1}^J \text{Normal}(\theta_j \mid \hat{\theta}_j, V_{\theta_j})$
Evaluated at $\theta_j = \hat{\theta}_j$, we get the marginal:
$p(\mu, \log \sigma, \log \tau \mid y) \propto \tau \left[ \prod_{j=1}^J \text{Normal}\big(\hat{\theta}_j \mid \mu, \tau^2\big) \right] \left[ \prod_{j=1}^J \prod_{i=1}^{n_j} \text{Normal}\big(y_{ij} \mid \hat{\theta}_j, \sigma^2\big) \right] \left[ \prod_{j=1}^J V_{\theta_j}^{1/2} \right]$
The last product $\prod_j V_{\theta_j}^{1/2}$ comes from the normalizing constants of the normal densities in the denominator (each contributes a factor proportional to $V_{\theta_j}^{-1/2}$, and taking the ratio inverts that).
This is an exact expression for the marginal posterior of $(\mu, \log \sigma, \log \tau)$, up to a multiplicative constant independent of these parameters.
6. Finding the marginal posterior mode of (μ, log σ, log τ) via EM
Directly maximizing is messy because:
- The $\hat{\theta}_j$’s and $V_{\theta_j}$’s depend on $\mu, \sigma, \tau$,
- So the expression has a complicated dependence on $(\mu, \sigma, \tau)$.
Instead, we use the EM algorithm, treating:
- $\phi = (\mu, \log \sigma, \log \tau)$ as “parameters of interest”,
- $\gamma = \theta$ as “latent variables” (or “missing data”).
6.1 Joint log posterior
The joint log posterior (up to a constant) is:
$\begin{aligned} \log p(\theta, \mu, \log \sigma, \log \tau \mid y) &= -n \log \sigma – (J – 1)\log \tau \\ &\quad – \frac{1}{2\tau^2} \sum_{j=1}^J (\theta_j – \mu)^2 – \frac{1}{2\sigma^2} \sum_{j=1}^J \sum_{i=1}^{n_j} (y_{ij} – \theta_j)^2 + \text{const}, \end{aligned}$ (13.13)
where $n = \sum_j n_j$ is the total number of observations.
6.2 E-step: expectations over θ
In the E-step, we take the expectation of this log posterior with respect to the conditional posterior of θ given the current parameter values $(\mu^{\text{old}}, \sigma^{\text{old}}, \tau^{\text{old}})$.
We need:
- $E_{\text{old}}\big[(\theta_j – \mu)^2\big]$
- $E_{\text{old}}\big[(y_{ij} – \theta_j)^2\big]$
Using the conditional posterior:
$\theta_j \mid \mu, \sigma, \tau, y \sim \text{Normal}\big(\hat{\theta}_j, V_{\theta_j}\big)$
we get:
- For each j,
- $E_{\text{old}}\big[(\theta_j – \mu)^2\big] = \big(E_{\text{old}}(\theta_j – \mu)\big)^2 + \text{Var}_{\text{old}}(\theta_j) = (\hat{\theta}_j – \mu)^2 + V_{\theta_j}$
- where $\hat{\theta}_j, V_{\theta_j}$ are computed using the old values $(\mu^{\text{old}}, \sigma^{\text{old}}, \tau^{\text{old}})$.
- Similarly, $E_{\text{old}}\big[(y_{ij} – \theta_j)^2\big] = (y_{ij} – \hat{\theta}_j)^2 + V_{\theta_j}.$
Plugging these expectations back into (13.13) gives the expected log posterior,
$E_{\text{old}}[\log p(\theta, \mu, \log \sigma, \log \tau \mid y)]$, as a function of $(\mu, \sigma, \tau)$ only.
6.3 M-step: maximizing over (μ, σ, τ)
In the M-step, we maximize that expected log posterior with respect to $\mu, \sigma, \tau$. The maximizing values are:
$\begin{aligned} \mu^{\text{new}} &= \frac{1}{J} \sum_{j=1}^J \hat{\theta}_j, \\[4pt] \sigma^{\text{new}} &= \left\{\frac{1}{n} \sum_{j=1}^J \sum_{i=1}^{n_j} \big[(y_{ij} – \hat{\theta}_j)^2 + V_{\theta_j}\big]\right\}^{1/2}, \\[4pt] \tau^{\text{new}} &= \left\{\frac{1}{J – 1} \sum_{j=1}^J \big[(\hat{\theta}_j – \mu^{\text{new}})^2 + V_{\theta_j}\big]\right\}^{1/2}. \end{aligned}$ (13.14)
These come from setting derivatives of the E-step objective to zero (details are in the exercises in the book).
We repeat E- and M-steps until convergence. Each EM iteration increases $\log p(\mu, \log \sigma, \log \tau \mid y)$, so we converge to a local marginal mode of the 3-dimensional marginal posterior.
In the coagulation example:
- Starting from the joint mode as an initial point, the EM algorithm converges in about 3 iterations.
- As is typical in hierarchical models, the marginal mode has larger estimates of σ and τ than the joint mode, because integrating over the $\theta_j$’s increases uncertainty.
7. Constructing a normal approximation to the joint posterior
Once we have:
- The marginal mode $(\hat{\mu}, \log \hat{\sigma}, \log \hat{\tau})$,
- And the 3×3 Hessian (second derivative) matrix of $\log p(\mu, \log \sigma, \log \tau \mid y)$ evaluated at this mode,
we can build a multivariate normal approximation:
$(\mu, \log \sigma, \log \tau) \mid y \approx \mathcal{N}\big( (\hat{\mu}, \log \hat{\sigma}, \log \hat{\tau}), V_\phi \big)$
where:
- $V_\phi = \left(-\nabla^2 \log p(\mu, \log \sigma, \log \tau \mid y)\big|_{\text{mode}}\right)^{-1}$
To simulate from an approximation to the full joint posterior:
- Draw $(\mu^{(s)}, \log \sigma^{(s)}, \log \tau^{(s)})$ from this multivariate normal.
- Propose $\sigma^{(s)} = \exp(\log \sigma^{(s)})$, $\tau^{(s)} = \exp(\log \tau^{(s)})$.
- For each draw, sample:
- $\theta_j^{(s)} \sim \text{Normal}\big(\hat{\theta}_j(\mu^{(s)}, \sigma^{(s)}, \tau^{(s)}), \, V_{\theta_j}(\mu^{(s)}, \sigma^{(s)}, \tau^{(s)})\big)$,
- using the exact conditional posterior $p(\theta \mid \mu, \sigma, \tau, y)$.
This gives a set of approximate draws from $p(\theta, \mu, \sigma, \tau \mid y)$.
The Table reports posterior quantiles (2.5%, 25%, median, 75%, 97.5%) for:
- Each $θ_j$
- μ, σ, τ
based on this normal approximation in $(\mu, \log \sigma, \log \tau)$ and exact conditional sampling for θ.
8. Using the approximation in practice
If the resulting approximate posterior looks reasonable (e.g., compared with MCMC or prior expectations), we can:
- Use it directly for posterior means, medians, and intervals.
- Use it as a starting point or proposal distribution for more advanced algorithms (e.g., importance sampling, MCMC, variational methods).
If the approximation is not adequate (e.g., due to strong skewness or multimodality):
- The normal approximation is still useful as:
- A baseline comparison for more advanced computation (like Gibbs, HMC, Stan),
- A provider of good initial values or over-dispersed proposals (e.g., switching to a multivariate t with low degrees of freedom and then using importance resampling).