Consider a dataset of size $n$ with $d$-dimensional multivariate observations drawn from
$N_d(\mu, \Sigma)$, where each observation is partitioned into observed components $y_{\text{obs}}$ and missing components $y_{\text{mis}}$. Methods are presented for a uniform prior on $(\mu, \Sigma)$, with extensions to hierarchical models introduced in the next section.

1. Posterior Mode via EM

The multivariate normal is an exponential family with sufficient statistics

  • $\sum_{i=1}^n y_{ij}$ for $j = 1,\dots,d$
  • $\sum_{i=1}^n y_{ij} y_{ik}$ for $j,k = 1,\dots,d$

Let $y^{\text{obs}}_i$ denote observed components of $y_i$ and $y^{\text{mis}}_i$ the missing components.
Let $\theta^{\text{old}} = (\mu^{\text{old}}, \Sigma^{\text{old}})$ denote current parameter estimates.

E-step

Compute conditional expectations of sufficient statistics:∑i=1nE(yij∣yobs,θold)=∑i=1nyijold,\sum_{i=1}^n \mathbb{E}(y_{ij} \mid y_{\text{obs}}, \theta^{\text{old}}) = \sum_{i=1}^n y^{\text{old}}_{ij},∑i=1nE(yijyik∣yobs,θold)=∑i=1n(yijoldyikold+cijkold),\sum_{i=1}^n \mathbb{E}(y_{ij} y_{ik} \mid y_{\text{obs}}, \theta^{\text{old}}) = \sum_{i=1}^n \left(y^{\text{old}}_{ij} y^{\text{old}}_{ik} + c^{\text{old}}_{ijk}\right),

where:yijold={yij,yij observed,E(yij∣yobs,θold),yij missing,y^{\text{old}}_{ij}= \begin{cases} y_{ij}, & y_{ij} \text{ observed},\\[4pt] \mathbb{E}(y_{ij}\mid y_{\text{obs}}, \theta^{\text{old}}), & y_{ij} \text{ missing}, \end{cases}cijkold={0,yij or yik observed,cov⁡(yij,yik∣yobs,θold),yij,yik missing.c^{\text{old}}_{ijk} = \begin{cases} 0, & y_{ij} \text{ or } y_{ik} \text{ observed},\\[4pt] \operatorname{cov}(y_{ij}, y_{ik}\mid y_{\text{obs}}, \theta^{\text{old}}), & y_{ij}, y_{ik} \text{ missing}. \end{cases}

The conditional mean vector and covariance matrix for $y^{\text{mis}}_i$ follow from multivariate normal conditioning formulas.

M-step

Update parameters using expected complete-data statistics:μjnew=1n∑i=1nyijold,\mu^{\text{new}}_j = \frac{1}{n} \sum_{i=1}^n y^{\text{old}}_{ij},σjknew=1n∑i=1n(yijoldyikold+cijkold−μjnewμknew).\sigma^{\text{new}}_{jk} = \frac{1}{n} \sum_{i=1}^n \left( y^{\text{old}}_{ij} y^{\text{old}}_{ik} + c^{\text{old}}_{ijk} – \mu^{\text{new}}_j \mu^{\text{new}}_k \right).

The initial estimate of $\Sigma$ must be positive definite; complete-case estimates are often useful.

2. Sampling from the Posterior Distribution

Starting from EM mode estimates, posterior draws for $(y_{\text{mis}}, \mu, \Sigma)$ can be obtained using data augmentation:

  1. Draw $y_{\text{mis}}$ given $(\mu, \Sigma)$.
  2. Draw $\mu$ and $\Sigma$ given completed data.

This reproduces the standard Gibbs sampler for the multivariate normal.
If needed, Metropolis updates replace steps not available in closed form.

Efficiency under monotone missingness

If variables can be ordered so that observed-patterns follow a monotone structure, the posterior for $\theta$ can be sampled directly using a blockwise regression-factorization. Let $\psi = \psi(\theta) = (\psi_1, \dots, \psi_k)$, where each block corresponds to conditional regression parameters for a segment in the monotone pattern.

The log-likelihood factorizes:log⁡p(yobs∣ψ)=log⁡p(y1obs∣ψ1)+log⁡p(y2obs∣y1obs,ψ2)+⋯+log⁡p(ykobs∣y1:k−1obs,ψk).\log p(y_{\text{obs}} \mid \psi) = \log p(y^{\text{obs}}_1 \mid \psi_1) + \log p(y^{\text{obs}}_2 \mid y^{\text{obs}}_1, \psi_2) + \cdots + \log p(y^{\text{obs}}_k \mid y^{\text{obs}}_{1:k-1}, \psi_k).

If the prior similarly factorizes, posterior draws can be obtained sequentially:ψ1∣yobs,ψ2∣ψ1,yobs,…ψk∣ψ1:(k−1),yobs.\psi_1 \mid y_{\text{obs}}, \quad \psi_2 \mid \psi_1, y_{\text{obs}}, \quad \dots \quad \psi_k \mid \psi_{1:(k-1)}, y_{\text{obs}}.

For nearly monotone patterns, a partial imputation step can be added to create a monotone structure, improving efficiency over full data augmentation.

3. Extending to the Multivariate t Model

Following the scale-mixture representation, assume:yi∣θ,Vi∼Nd(μ,ViΣ),y_i \mid \theta, V_i \sim N_d(\mu, V_i \Sigma),

with latent scale parametersVi∼Inv-χ2(ν,1),V_i \sim \text{Inv-}\chi^2(\nu, 1),

where $\nu$ is fixed.

Data augmentation

Each iteration includes:

  1. Imputing $V_i$

p(Vi∣yiobs,θ,ν)∝N ⁣(yiobs∣μiobs,ViΣiobs)⋅Inv-χ2(Vi∣ν,1),p(V_i \mid y^{\text{obs}}_i, \theta, \nu) \propto N\!\left(y^{\text{obs}}_i \mid \mu^{\text{obs}}_i, V_i \Sigma^{\text{obs}}_i\right) \cdot \text{Inv-}\chi^2(V_i \mid \nu, 1),

which remains scaled inverse-$\chi^2$, so draws are straightforward.

  1. Imputing $y^{\text{mis}}_i$

Given $(V_i, \theta)$, sampling $y^{\text{mis}}_i$ follows the conditional multivariate normal distribution.

  1. Posterior update

Treating $(y^{\text{mis}}, V)$ as observed, update $(\mu, \Sigma)$ via weighted complete-data formulas.

EM version

Replace draws with expectations:

  • $\mathbb{E}(V_i)$ from the inverse-$\chi^2$ conditional
  • conditional means and covariances for $y^{\text{mis}}_i$

The M-step maximizes the posterior rather than sampling.

If $\nu$ is also unknown, ECME iterations update $\nu$ using the methods described in earlier chapters.

4. Nonignorable Missing-Data Models

For nonignorable mechanisms, inference conditions on both $(y_{\text{obs}}, I)$ and includes parameters $\phi$ for the missingness process. EM or Gibbs sampling proceeds as in the ignorable case, but each step explicitly includes the missing-data mechanism:

  • parameter updates involve $(\theta, \phi)$
  • conditional distributions depend on $(y_{\text{obs}}, I)$
  • expectations or draws for $y_{\text{mis}}$ incorporate the nonignorable mechanism