Conjugate Analysis for the Multivariate Gaussian

The goal of this post is to derive the likelihood, posterior, and posterior predictive for a multivariate Gaussian model with an unknown mean parameter. Consider the $D$-variate Gaussian,

$$\mathbf{x}\sim {\mathcal{N}}_D(\mathbf{\mu },\mathbf{\Sigma }),\text{\hspace{1em}(1)}$$

with density function $f(\cdot )$:

$$f(\mathbf{x}\mid \mathbf{\mu },\mathbf{\Sigma })=(2\pi {)}^{-D/2}\det (\mathbf{\Sigma }{)}^{-1/2}\exp \left(-\frac{1}{2}(\mathbf{x}-\mathbf{\mu }{)}^{\top }{\mathbf{\Sigma }}^{-1}(\mathbf{x}-\mathbf{\mu })\right).\text{\hspace{1em}(2)}$$

The likelihood is over $N$ i.i.d. observations, denoted with the design matrix $\mathbf{X}$, is

$$\begin{array}{rl}\mathcal{L}(\mathbf{X}\mid \mathbf{\mu },\mathbf{\Sigma })& =\prod _{n=1}^{N}f({\mathbf{x}}_n\mid \mathbf{\mu },\mathbf{\Sigma })\\ & =(2\pi {)}^{-ND/2}\det (\mathbf{\Sigma }{)}^{-N/2}\exp \left(-\frac{1}{2}\sum _{n=1}^N({\mathbf{x}}_n-\mathbf{\mu }{)}^{\top }{\mathbf{\Sigma }}^{-1}({\mathbf{x}}_n-\mathbf{\mu })\right).\end{array}\text{\hspace{1em}(3)}$$

Assume $\mathbf{\Sigma }$ is known. A common prior for the mean parameter $\mathbf{\mu }$ is another Gaussian,

$$\pi (\mathbf{\mu })={\mathcal{N}}_D({\mathbf{\mu }}_0,{\mathbf{\Sigma }}_0).\text{\hspace{1em}(4)}$$

Our goal is to show that $\pi (\mathbf{\mu })$ is a conjugate prior, meaning the posterior $p(\mathbf{\mu }\mid \mathbf{x})$ is also Gaussian. The posterior is

$$p(\mathbf{\mu }\mid \mathbf{X})\propto \mathcal{L}(\mathbf{X}\mid \mathbf{\mu },\mathbf{\Sigma })\pi (\mathbf{\mu }),\text{\hspace{1em}(5)}$$

which we can write out explicitly as

$$\begin{array}{rl}& \mathcal{L}(\mathbf{X}\mid \mathbf{\mu },\mathbf{\Sigma })\pi (\mathbf{\mu })\\ & \propto \exp \left(-\frac{1}{2}\sum _{n=1}^N({\mathbf{x}}_n-\mathbf{\mu }{)}^{\top }{\mathbf{\Sigma }}^{-1}({\mathbf{x}}_n-\mathbf{\mu })\right)\exp \left(-\frac{1}{2}(\mathbf{\mu }-{\mathbf{\mu }}_0{)}^{\top }{\mathbf{\Sigma }}_0^{-1}(\mathbf{\mu }-{\mathbf{\mu }}_0)\right).\end{array}\text{\hspace{1em}(6)}$$

We can write the terms inside the exponents as

$$\sum _{n=1}^N{\mathbf{x}}_n^{\top }{\mathbf{\Sigma }}^{-1}{\mathbf{x}}_n+N{\mathbf{\mu }}^{\top }{\mathbf{\Sigma }}^{-1}\mathbf{\mu }-2N\bar{\mathbf{x}}{\mathbf{\Sigma }}^{-1}\mathbf{\mu }+{\mathbf{\mu }}^{\top }{\mathbf{\Sigma }}_0^{-1}\mathbf{\mu }+{\mathbf{\mu }}_0^{\top }{\mathbf{\Sigma }}_0^{-1}{\mathbf{\mu }}_0-2{\mathbf{\mu }}^{\top }{\mathbf{\Sigma }}_0^{-1}{\mathbf{\mu }}_0,\text{\hspace{1em}(7)}$$

where $\bar{\mathbf{x}}=1/N{\sum }_{n=1}^N{\mathbf{x}}_n$. Since the posterior is w.r.t. to $\mathbf{\mu }$, we can drop terms that do not depend on $\mathbf{\mu }$—we’ll still retain the Gaussian kernel and can just properly normalize it—and combine like terms:

$$p(\mathbf{\mu }\mid \mathbf{X})\propto \exp \left(-\frac{1}{2}\left[{\mathbf{\mu }}^{\top }(N{\mathbf{\Sigma }}^{-1}+{\mathbf{\Sigma }}_0^{-1})\mathbf{\mu }-2\mathbf{\mu }(N{\mathbf{\Sigma }}^{-1}\bar{\mathbf{x}}+{\mathbf{\Sigma }}_0^{-1}{\mathbf{\mu }}_0)\right]\right).\text{\hspace{1em}(8)}$$

Now we just complete the square. First, let

$$\begin{array}{rl}\mathbf{M}& =N{\mathbf{\Sigma }}^{-1}+{\mathbf{\Sigma }}_0^{-1},\\ \mathbf{b}& =N{\mathbf{\Sigma }}^{-1}\bar{\mathbf{x}}+{\mathbf{\Sigma }}_0^{-1}{\mathbf{\mu }}_0.\end{array}\text{\hspace{1em}(9)}$$

Then Eq. $8$ is equivalent to

$$\begin{array}{rl}p(\mathbf{\mu }\mid \mathbf{X})& \propto \exp \left(-\frac{1}{2}\left[(\mathbf{\mu }-{\mathbf{M}}^{-1}\mathbf{b}{)}^{\top }\mathbf{M}(\mathbf{\mu }-{\mathbf{M}}^{-1}\mathbf{b})-{\mathbf{b}}^{\top }{\mathbf{M}}^{-1}\mathbf{b}\right]\right),\\ & \Downarrow \\ p(\mathbf{\mu }\mid \mathbf{X})& ={\mathcal{N}}_D(\mathbf{\mu }\mid {\mathbf{M}}^{-1}\mathbf{b},{\mathbf{M}}^{-1}).\end{array}\text{\hspace{1em}(10)}$$

This is the posterior for a multivariate Gaussian with unknown mean. We can write it a standard form, e.g. see (Murphy, 2007), as:

$$\begin{array}{rl}p(\mathbf{\mu }\mid \mathbf{X})& ={\mathcal{N}}_D(\mathbf{\mu }\mid {\mathbf{\mu }}_N,{\mathbf{\Sigma }}_N),\\ {\mathbf{\Sigma }}_N& =(N{\mathbf{\Sigma }}^{-1}+{\mathbf{\Sigma }}_0^{-1}{)}^{-1},\\ {\mathbf{\mu }}_N& ={\mathbf{\Sigma }}_N(N{\mathbf{\Sigma }}^{-1}\bar{\mathbf{x}}+{\mathbf{\Sigma }}_0^{-1}{\mathbf{\mu }}_0).\end{array}\text{\hspace{1em}(11)}$$

To compute the posterior predictive,

$$p({\mathbf{x}}_{\text{new}}\mid \mathbf{X})=\int f(\mathbf{x}\mid \mathbf{\mu },\mathbf{\Sigma })p(\mathbf{\mu }\mid \mathbf{X})d\mathbf{\mu },\text{\hspace{1em}(12)}$$

we observe the following:

$$\begin{array}{rl}({\mathbf{x}}_{\text{new}}-\mathbf{\mu })& \sim \mathcal{N}({\mathbf{x}}_{\text{new}}-\mathbf{\mu }\mid 0,\mathbf{\Sigma }),\\ \mathbf{\mu }& \sim \mathcal{N}(\mathbf{\mu }\mid {\mathbf{\mu }}_N,{\mathbf{\Sigma }}_N).\end{array}\text{\hspace{1em}(13)}$$

Notice that $({\mathbf{x}}_{\text{new}}-\mathbf{\mu })$ and $\mathbf{\mu }$ are independent. Intuitively, if I tell you the value of $\mathbf{\mu }$, that tells you nothing about the value of $({\mathbf{x}}_{\text{new}}-\mathbf{\mu })$ because the distribution on $({\mathbf{x}}_{\text{new}}-\mathbf{\mu })$ does not contain the parameter $\mathbf{\mu }$. Formally,

$$p({\mathbf{x}}_{\text{new}}-\mathbf{\mu }\mid \mathbf{\mu })=p({\mathbf{x}}_{\text{new}}-\mathbf{\mu }).\text{\hspace{1em}(14)}$$

Since $({\mathbf{x}}_{\text{new}}-\mathbf{\mu })$ and $\mathbf{\mu }$ are independent, we can simply add the means and covariances of the two distributions when adding the random variables, implying:

$$\begin{array}{rl}{\mathbf{x}}_{\text{new}}& =({\mathbf{x}}_{\text{new}}-\mathbf{\mu })+\mathbf{\mu }\\ & \Downarrow \\ p({\mathbf{x}}_{\text{new}}\mid \mathbf{X})& =\mathcal{N}({\mathbf{x}}_{\text{new}}\mid {\mathbf{\mu }}_N,\mathbf{\Sigma }+{\mathbf{\Sigma }}_N).\end{array}\text{\hspace{1em}(15)}$$