1. Goal: marginal posterior modes for a subset of parameters

We split the full parameter vector into two parts:

  • Full parameter vector: $\theta = (\gamma, \phi)$
  • $\phi$: the subset we care about (e.g., hyperparameters, group-level parameters)
  • $\gamma$: the remaining part, which we will treat like “latent / missing” variables

Often:

  • A joint normal approximation for all of θ is not useful (high dimension, non-normal joint shape).
  • But a good approximation to the marginal posterior of φ, $p(\phi \mid y)$, is useful.

Strategy:

  1. Use EM (or its variants) to find a marginal posterior mode of $p(\phi \mid y)$.
  2. Approximate $p(\phi \mid y)$ around that mode with:
    • a normal distribution, or
    • a t-distribution, or
    • a mixture of normals/t’s.
  3. Then approximate conditional $p(\gamma \mid \phi, y)$ (next section in the book), often as normal / t / mixture with parameters depending on φ.

In this section, they focus on step 1: finding marginal posterior modes of φ using EM, thinking of:

  • $\phi$ = “parameters we want to optimize over”
  • $\gamma$ = “missing data / auxiliary variables”

2. EM as a mode-finding algorithm for $p(\phi \mid y)$

Interpretation:

  • Treat $\phi$ as “parameters” and $\gamma$ as “missing data”.
  • The classic idea: alternate between:
    1. E-step: fill in missing information (take expectations of functions of γ),
    2. M-step: re-estimate φ as if those expected sufficient statistics were observed.

EM formalizes this as an algorithm that increases the marginal log posterior $\log p(\phi \mid y)$ at each iteration, and converges to a local mode of $p(\phi \mid y)$.

3. Deriving the EM / GEM update

We want to maximize the marginal posterior:

$p(\phi \mid y) = \int p(\gamma, \phi \mid y)\, d\gamma$

Basic identity:

$\log p(\phi \mid y) = \log p(\gamma, \phi \mid y) – \log p(\gamma \mid \phi, y)$

Now take expectations over γ using the current conditional distribution:

  • Let $\phi_{\text{old}}$​ be the current value,
  • Let $p(\gamma \mid \phi_{\text{old}}, y)$ be the current conditional posterior of γ,
  • Define: $E_{\text{old}}[\cdot] = \mathbb{E}[\cdot \mid \phi_{\text{old}}, y]$

Take expectation of both sides:

$\log p(\phi \mid y) = E_{\text{old}}[\log p(\gamma, \phi \mid y)] – E_{\text{old}}[\log p(\gamma \mid \phi, y)]$ (13.5)

Important facts:

  1. The left-hand side $\log p(\phi \mid y)$ does not depend on γ, so expectation over γ does nothing to it.
  2. The second term, $E_{\text{old}}[\log p(\gamma \mid \phi, y)],$ is maximized at $\phi = \phi_{\text{old}}$​. (They refer to an exercise for this.)

So if we choose $\phi_{\text{new}}$​ such that:

$E_{\text{old}}[\log p(\gamma, \phi_{\text{new}} \mid y)] > E_{\text{old}}[\log p(\gamma, \phi_{\text{old}} \mid y)]$

then:

  • The first term on the RHS of (13.5) increases,
  • The second term does not increase (it’s maximal at $\phi_{\text{old}}$​),

Therefore:

$\log p(\phi_{\text{new}} \mid y) > \log p(\phi_{\text{old}} \mid y)$

This is the key insight:

if you increase that expected complete-data log posterior, you increase the marginal log posterior.

Q-function and EM / GEM

Define the Q-function:

$Q(\phi \mid \phi_{\text{old}}) = E_{\text{old}}[\log p(\gamma, \phi \mid y)] = \int \log p(\gamma, \phi \mid y)\, p(\gamma \mid \phi_{\text{old}}, y)\, d\gamma$

Then:

  • EM algorithm: at each iteration, choose $\phi_{\text{new}}$​ that maximizes $Q(\phi \mid \phi_{\text{old}})$.
  • GEM (Generalized EM): choose $\phi_{\text{new}}$​ that merely increases $Q(\phi \mid \phi_{\text{old}})$ (not necessarily the maximum).

Both guarantee:

  • $\log p(\phi \mid y)$ increases every iteration,
  • EM therefore converges to a local mode of $p(\phi \mid y)$ (except in special edge cases),
  • GEM increases $\log p(\phi \mid y)$ but might not fully converge to a mode if you don’t maximize enough.

The speed of convergence depends on how much of the information about φ is “missing” when you integrate over γ. If a lot of information is missing (relative to having complete γ), EM can converge slowly.

4. EM algorithm in practice (φ–γ view)

Algorithm:

  1. Initialization
    Choose a crude starting value $\phi_0$​.
  2. For $t = 1, 2, \dots$:
    • (a) E-step
      • Compute the expected log joint posterior of $(\gamma, \phi)$, treating γ as random and averaging over $p(\gamma \mid \phi_{\text{old}}, y)$, where $\phi_{\text{old}} = \phi_{t-1}$​:
        • $Q(\phi \mid \phi_{\text{old}}) = E_{\text{old}}[\log p(\gamma, \phi \mid y)] = \int \log p(\gamma, \phi \mid y)\, p(\gamma \mid \phi_{\text{old}}, y)\, d\gamma$
    • (b) M-step
      • Choose $\phi_t$​ to maximize $Q(\phi \mid \phi_{\text{old}})$ (or at least increase it for GEM).

This guarantees:

  • $p(\phi_t \mid y)$ increases,
  • EM converges to a local marginal mode of $\phi$.

Multiple modes & debugging

  • To search for multiple modes, run EM from many starting points $\phi_0$​ spread out in parameter space.
  • After finding several modes, use a normal approximation around each mode to estimate their relative importance (as in the previous section on mixture approximations).
  • Debugging tip: compute $\log p(\phi_t \mid y)$ at each iteration and verify it is monotone increasing.

5. Example: Normal model with unknown mean and variance

Model:

  • Observations: $y_1, \dots, y_n$​ i.i.d.:
    • $y_i \sim N(\mu, \sigma^2)$
  • Prior:
    • $\mu \sim N(\mu_0, \tau_0^2), \quad \log \sigma \sim \text{Uniform} \quad (\text{noninformative on log scale})$

Goal:

  • Find the marginal posterior mode of μ, averaging over σ.
  • So here:
    • $\phi = \mu$ (parameter of interest),
    • $\gamma = \sigma$ (treated as “missing/latent” from EM’s point of view).

Joint log posterior

Ignoring constants not depending on μ or σ²:

$\log p(\mu, \sigma \mid y) = – \frac{1}{2\tau_0^2}(\mu – \mu_0)^2 – (n+1)\log \sigma – \frac{1}{2\sigma^2} \sum_{i=1}^n (y_i – \mu)^2 + \text{constant}$ (13.6)

E-step: expectation over σ given $μ_{\text{old}}$

We need:

$E_{\text{old}}[\log p(\mu, \sigma \mid y)]$

where the expectation is taken over $\sigma$ given $\mu_{\text{old}}$ and $y$.

If we plug (13.6) into this expectation, we get:

$E_{\text{old}}[\log p(\mu, \sigma \mid y)] = – \frac{1}{2\tau_0^2}(\mu – \mu_0)^2 – \frac{1}{2} E_{\text{old}}\left(\frac{1}{\sigma^2}\right) \sum_{i=1}^n (y_i – \mu)^2 – (n+1) E_{\text{old}}[\log \sigma] + \text{constant}$ (13.7)

The term $E_{\text{old}}[\log \sigma]$ does not depend on μ, so it doesn’t affect the M-step (when we optimize over μ). We only really need:

$E_{\text{old}}\left(\frac{1}{\sigma^2}\right)$

Conditional posterior for $\sigma^2$ given μ and y:

  • For fixed μ, this is the usual normal model with known mean and unknown variance, so:
    • $\sigma^2 \mid \mu, y \sim \text{Inv-}\chi^2\!\left(n,\; \frac{1}{n} \sum_{i=1}^n (y_i – \mu)^2 \right)$ (scaled inverse-chi-squared).

Thus:

  • $1/\sigma^2$ follows a scaled $\chi^2$,
  • Its expectation, given $\mu_{\text{old}}$​, is:
    • $E_{\text{old}}\left(\frac{1}{\sigma^2}\right) = \frac{n}{\sum_{i=1}^n (y_i – \mu_{\text{old}})^2}$

So (13.7) becomes:

$E_{\text{old}}[\log p(\mu, \sigma \mid y)] = – \frac{1}{2\tau_0^2}(\mu – \mu_0)^2 – \frac{1}{2} \left(\frac{n}{\sum_{i=1}^n (y_i – \mu_{\text{old}})^2}\right) \sum_{i=1}^n (y_i – \mu)^2 + \text{const}$ (13.8)

Interpretation of E-step result

Equation (13.8) looks like the log posterior of a normal model:

  • Prior: $\mu \sim N(\mu_0, \tau_0^2)$
  • “Data”: $y_i$​ each with effective variance:
    • $\frac{1}{E_{\text{old}}(1/\sigma^2)} = \frac{\sum_{i=1}^n (y_i – \mu_{\text{old}})^2}{n}$

So in the M-step, μ gets updated like a posterior mean/mode of a normal-normal model with this variance.

M-step: updating μ

Maximizing (13.8) with respect to μ gives the M-step update:

$\mu_{\text{new}} = \frac{ \frac{1}{\tau_0^2}\mu_0 + \left(\frac{n}{\sum_{i=1}^n (y_i – \mu_{\text{old}})^2}\right) \sum_{i=1}^n y_i }{ \frac{1}{\tau_0^2} + \frac{n^2}{\sum_{i=1}^n (y_i – \mu_{\text{old}})^2} }$

(Their expression is equivalent, written slightly differently.)

Iterating this:

  • Start from some $\mu_0$​,
  • Repeat E-step and M-step,
  • $\mu_t$​ converges to the marginal posterior mode of $p(\mu \mid y)$, integrating out σ.

6. Extensions of EM: ECM, ECME, AECM

The basic EM can be slow. Several extensions generalize the M-step to speed up convergence or make computations easier.

6.1 ECM (Expectation–Conditional Maximization)

Idea:

  • Instead of maximizing over the full φ at once in the M-step, we:
    • Partition φ into S blocks: $\phi = (\phi_1, \dots, \phi_S)$
    • Perform conditional maximizations one block at a time.

Algorithm (one EM iteration):

  1. E-step: compute $Q(\phi \mid \phi_{\text{old}})$ as usual.
  2. CM-steps: for $s = 1, \dots, S$:
    • Hold all components except $\phi_s$​ fixed,
    • Maximize $Q(\phi \mid \phi_{\text{old}})$ with respect to $\phi_s$​ only.

This is still a GEM algorithm:

  • Each CM-step increases (or at least doesn’t decrease) the Q-function,
  • So the marginal posterior density $p(\phi \mid y)$ still increases each iteration,
  • If each CM-step is a true maximization (setting derivatives to zero), ECM also converges to a local mode.

Multicycle ECM: Do extra E-steps within one overall iteration (e.g., E–CM–E–CM–…); more work per iteration but sometimes fewer iterations.

6.2 ECME (Expectation / Conditional Maximization Either)

ECME extends ECM:

  • Some CM-steps maximize $Q(\phi \mid \phi_{\text{old}})$ (like ECM),
  • Other CM-steps maximize the actual log marginal posterior $\log p(\phi \mid y)$ directly.

Benefits:

  • Can use faster optimization methods (e.g., Newton’s method) on the actual marginal target for some blocks.
  • Often much faster convergence, because you are directly improving the real objective $\log p(\phi \mid y)$ in some steps.

Properties:

  • $\log p(\phi \mid y)$ still increases at each iteration,
  • With proper maximization (setting derivatives to zero), ECME converges to a local mode.

6.3 AECM (Alternating ECM)

AECM generalizes further:

  • Allow different definitions of γ (the “missing data”) at different CM-steps.
  • In other words, the augmented data structure can change between steps.
  • This makes it possible to choose the most convenient incomplete / complete data structure for each parameter block.

Benefits:

  • In some problems, this gives significantly faster convergence.
  • Still retains:
    • Monotone increase in the posterior density,
    • Convergence to a local mode.

7. Supplemented EM / ECM: approximating variance at the mode

EM (and ECM, ECME, AECM) give you:

  • A marginal posterior mode $\hat{\phi}$​,
  • But we also want an approximate covariance / variance matrix to build a normal approximation to $p(\phi \mid y)$.

The Supplemented EM (SEM) and Supplemented ECM (SECM) algorithms:

  • Use the EM (or ECM) mapping and derivatives to estimate the asymptotic variance at the mode.

7.1 SEM: notation and idea

Define:

  • EM mapping: $\phi_{t+1} = M(\phi_t)$
  • At the mode $\hat{\phi}$​, we care about the Jacobian of this map:
    • $D M = \left.\frac{\partial M(\phi)}{\partial \phi^T}\right|_{\phi = \hat{\phi}}$

Define:

  • $V_{\text{joint}}$​: asymptotic variance matrix based on the joint log posterior $\log p(\phi, \gamma \mid y)$ averaged over γ:
    • $V_{\text{joint}} = \left( – E \left[\frac{\partial^2 \log p(\phi, \gamma \mid y)}{\partial \phi\, \partial \phi^T} \right] \right)^{-1}\Bigg|_{\phi = \hat{\phi}}$

Then the marginal asymptotic variance of φ is:

$V = V_{\text{joint}} + V_{\text{joint}}\, D M \, (I – D M)^{-1}$

So SEM needs:

  • $V_{\text{joint}}$ (often can be derived analytically),
  • $D M$ (estimated numerically using the EM mapping near the mode).

7.2 Estimating DM numerically

At a high level:

  1. Run EM to convergence ⇒ get $\hat{\phi}$.
  2. Choose a starting φ₀ for SEM (not equal to $\hat{\phi}$​).
  3. For each iteration t (SEM iterations), compute a matrix $R_t$​ whose rows approximate rows of DMD MDM:
    • For each component i of φ:
      1. Construct a perturbed vector $\phi_t^{(i)}$​ equal to $\hat{\phi}$​ except the i-th element is replaced by $\phi_{t,i}$,
      2. Run one EM E-step+M-step from $\phi_t^{(i)}$ to get $\phi_{t+1}^{(i)}$​,
      3. Use finite differences between $\phi_{t+1}^{(i)}$​ and $\phi_t^{(i)}$ (relative to $\hat{\phi}$​) to estimate the row i of the Jacobian.

As $t$ grows, $R_t$​ converges to $D M$. Then you plug into the formula for V.

7.3 SECM

For ECM:

  • The convergence behavior is more complicated, because ECM uses multiple CM-steps.
  • The SECM formula uses:
    • $D^{\text{ECM}}_M$: Jacobian of the ECM mapping,
    • $D^{\text{CM}}_M$​: convergence matrix for conditional maximization applied directly to $\log p(\phi \mid y)$,
    • Constraint gradients $∇ g_s(\hat{\phi})$ for each CM-block.

The final formula:

$V = V_{\text{joint}} + V_{\text{joint}} (D^{\text{ECM}}_M – D^{\text{CM}}_M) (I – D^{\text{ECM}}_M)^{-1}$

where $D^{\text{CM}}_M$ is constructed from $V_{\text{joint}}$​ and the gradients of the constraint functions. This is more technical but the idea is similar: combine information from the EM mapping and the joint curvature.

8. Parameter-expanded EM (PX-EM)

Parameter expansion for EM is the mode-finding analogue of parameter-expanded Gibbs sampling.

Idea:

  • Introduce extra parameters (like α in the t-example of parameter expansion for Gibbs).
  • These parameters are not identified by the data and serve only to improve convergence.
  • In EM, we treat these extra parameters as part of φ, and use ECM to update them along with the original parameters.

Example (connection to earlier t-model):

  • In the t-model, we introduced latent variances $V_i$ and a scale parameter α, etc.
  • Originally, that parameter expansion was motivated in the EM context:
    • Treat the $V_i$ as missing data (γ),
    • Use ECM to maximize over $(\mu, \tau, \alpha)$ (φ),
    • The expansion reduces dependence and speeds up convergence of the EM iterations.

PX-EM and PX-Gibbs share the same logic:

Add non-identified parameters in a careful way to break strong dependencies and improve mixing (for Gibbs) or convergence (for EM) while ensuring that the marginal posterior / likelihood in the original parameter of interest is unchanged.