Multivariate Normal with Unknown Mean and Variance#

Part 1 · Stage 3 · 🧮 Multiparameter Models · Lesson 025 of 144 · beginner

◀ Previous · Multivariate Normal Model with Known Variance · Next · Example: Bayesian analysis of a bioassay experiment (logistic, nonconjugate) ▶ · ↑ 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 full multivariate model#

Drop the fiction that \(\Sigma\) is known. With \(y_i \sim \mathrm{N}(\theta, \Sigma)\) and both unknown, the conjugate prior is the matrix version of the normal–inverse-\(\chi^2\) from earlier in this stage: a normal–inverse-Wishart,

\[\Sigma \sim \text{Inv-Wishart}(\nu_0, \Lambda_0^{-1}), \qquad \theta \mid \Sigma \sim \mathrm{N}\!\left(\mu_0,\, \frac{\Sigma}{\kappa_0}\right).\]

Again the mean’s prior is scaled by the very covariance being estimated, and again \(\kappa_0\) and \(\nu_0\) are equivalent sample sizes — for the mean and for the covariance respectively.

The update#

Conjugacy holds, and the hyperparameters update by addition exactly as in the scalar case:

\[\kappa_n = \kappa_0 + n, \qquad \nu_n = \nu_0 + n, \qquad \mu_n = \frac{\kappa_0 \mu_0 + n \bar{y}}{\kappa_0 + n},\]

\[\Lambda_n = \Lambda_0 + S + \frac{\kappa_0 n}{\kappa_0 + n} (\bar{y} - \mu_0)(\bar{y} - \mu_0)^{\top},\]

where \(S = \sum_i (y_i - \bar{y})(y_i - \bar{y})^{\top}\). The last term is the multivariate conflict detector: prior–data disagreement about the mean inflates the posterior covariance. Marginalising \(\Sigma\) gives a multivariate \(t\) for \(\theta\), the vector analogue of the scalar result.

Why the inverse-Wishart disappoints#

The inverse-Wishart is conjugate, guarantees positive-definiteness, and makes Gibbs sampling trivial — which is why it dominated Bayesian software for decades. It is also, in Gelman’s own later assessment, a poor default. A single degrees-of-freedom parameter \(\nu_0\) controls the certainty of every variance at once; the implied marginal for each variance has little density near zero, biasing small variances upward; and it forces an a priori dependence between correlations and variances that nobody believes.

The modern alternative#

Modern practice decomposes \(\Sigma\) into scales and a correlation matrix, giving each its own prior — an LKJ prior on the correlations, half-normal or half-Cauchy on the standard deviations:

import pymc as pm
with pm.Model():
    sd = pm.HalfNormal("sd", sigma=1, shape=d)             # scales, separately
    chol, corr, _ = pm.LKJCholeskyCov("Sigma", n=d, eta=2,  # correlations, separately
                                       sd_dist=pm.HalfNormal.dist(1))
    pm.MvNormal("y", mu=theta, chol=chol, observed=Y)

Conjugacy bought tractability when computation was scarce. With HMC available, the honest prior wins — a theme that recurs whenever Part I’s closed forms meet Part III’s samplers.

Hint