# Estimating the covariance posterior distribution of a multivariate gaussian

I need to "learn" the distribution of a bivariate gaussian with few samples, but a good hypothesis on the prior distribution, so I would like to use the bayesian approach.

I defined my prior: $$\mathbf{P}(\mathbf{\mu}) \sim \mathcal{N}(\mathbf{\mu_0},\mathbf{\Sigma_0})$$ $$\mathbf{\mu_0} = \begin{bmatrix} 0 \\ 0 \end{bmatrix} \ \ \ \mathbf{\Sigma_0} = \begin{bmatrix} 16 & 0 \\ 0 & 27 \end{bmatrix}$$

And my distribution given the hypothesis $$\mathbf{P}(x|\mathbf{\mu},\mathbf{\Sigma}) \sim \mathcal{N}(\mathbf{\mu},\mathbf{\Sigma})$$ $$\mathbf{\mu} = \begin{bmatrix} 0 \\ 0 \end{bmatrix} \ \ \ \mathbf{\Sigma} = \begin{bmatrix} 18 & 0 \\ 0 & 18 \end{bmatrix}$$

Now I know thanks to here that to estimate the mean given the data

$$\mathbf{P} (\mathbf{\mu} | \mathbf{x_1}, \dots , \mathbf{x_n}) \sim \mathcal{N}(\mathbf{\hat{\mu}_n}, \mathbf{\hat{\Sigma}_n})$$

I can compute:

$$\mathbf{\hat{\mu}_n} = \mathbf{\Sigma_0} \left( \mathbf{\Sigma_0} + {1 \over n} \mathbf{\Sigma} \right) ^ {-1} \left( {1 \over n} \sum_{i=1}^{n} \mathbf{x_i} \right) + {1 \over n} \mathbf{\Sigma} \left( \mathbf{\Sigma_0} + {1 \over n} \mathbf{\Sigma} \right) ^{-1} \mathbf{\mu_0}$$

$$\mathbf {\hat{\Sigma}_n} = {1 \over n} \mathbf{\Sigma_0} \left( \mathbf{\Sigma_0} + {1 \over n} \mathbf{\Sigma} \right) ^{-1} \mathbf{\Sigma}$$

Now comes the question, maybe I'm wrong, but it seems to me that $\mathbf{\Sigma_n}$ is just the covariance matrix for the estimated parameter $\mathbf{\mu_n}$, and not the estimated covariance of my data. What I would like would be to compute also

$$\mathbf{P} (\mathbf{\Sigma_{n_1}} | \mathbf{x_1}, \dots , \mathbf{x_n})$$

in order to have a fully specified distribution learned from my data.

Is this possible? Is it already solved by computing $\mathbf{\Sigma_n}$ and it's just expressed in the wrong way the formula above (or I am simply misentrepreting it)? References would be appreciated. Thanks a lot.

EDIT

From the comments, it appeared that my approach was "wrong", in the sense that I was assuming a constant covariance, defined by $\mathbf{\Sigma}$. What I need would be to put a prior also on it, $\mathbf{P}(\mathbf{\Sigma})$, but I don't know what distribution I should use, and subsequently what is the procedure to update it.

-
You have already specified the covariance of your data as $\mathbf{\Sigma} = \begin{bmatrix} 18 & 0 \\ 0 & 18 \end{bmatrix}$ - and you haven't specified a prior distribution for that to be updated from? – Corone Feb 26 '13 at 8:28
I see your point. So with my approach I basically assumed that the variance was constant and specified. If I want to estimate it, I need a prior on it. $\mathbf{P}(\mathbf{\Sigma}) \sim \mathcal{F} (\mathbf{\mu_{\Sigma}} , \Sigma_{\Sigma})$ Now, my problem is that it's not clear how to define it, and what would be an appropriate distribution for it, but this seems to be out of the scope of the first question. – unziberla Feb 26 '13 at 14:12
Then change the question :-) – Corone Feb 26 '13 at 14:20

You can do Bayesian updating for the covariance structure in much the same spirit as you updated the mean. The conjugate prior for the covariance matrix of the multivariate-normal is the Inverse-Wishart distribution, so it makes sense to start there,

$P(\Sigma) \sim W^{-1}(\mathbf{\Psi}, \nu)$

Then when you get your sample $X$ of length $n$ you can calculate the sample covariance estimate $\Sigma_X = \frac{1}{n}(X-\mu)^\top(X-\mu)$

This can then be used to update your estimate of the covariance matrix

$P(\Sigma|X) \sim W^{-1}(n\Sigma_X + \mathbf{\Psi}, n + \nu)$

You may choose to use the mean of this as your point estimate for the covariance (Posterior Mean Estimator)

$E[\Sigma|X] = \frac{n\Sigma_X + \mathbf{\Psi}}{\nu+n-p-1}$

or you might choose to use the mode (Maximum A Posteriori Estimator)

$\text{Mode}[\Sigma|X] = \frac{n\Sigma_X + \mathbf{\Psi}}{\nu+n+p+1}$

-
Thanks a lot. Now I assume something will change in my estimation process. As a first step, I should estimate the covariance $\mathbf{\hat{\Sigma}}$ with your procedure, then my distribution given the estimated hypothesis woulb be $\mathbf{P} (\mathbf{X} | \mu, \mathbf{\hat{\Sigma}} )$ and since $\mathbf{\hat{\Sigma}}$ is estimated and has its own distribution I am pretty sure this will somehow change my previous formula to compute $\mathbf{\hat{\mu}_n}$ (as it happens on gaussian MLE when using the sample variance). – unziberla Feb 26 '13 at 18:04
The approach that you describe would be instead to use $\mathbf{\hat{\Sigma}} = E[ \Sigma | \mathbf{x_1} \dots \mathbf{x_n} ]$ so that I have an actual value for the covariance, as if I knew it before. In a frequentist approach, this would sound wrong, but maybe there is something that I am missing from the fact that I assume the prior is known and this makes the procedure correct? – unziberla Feb 26 '13 at 18:05

Ok, I found the real solution for my problem. I am posting it even if the correct answer to my (misplaced) question is the one selected.

Basically, my question explains how to estimate the mean knowing the covariance, and the answer how to estimate the covariance knowing the mean. But my actual problem was estimating with both parameters unknown.

I found the answer on Wikipedia with the derivation explained here. The multivariate normal's conjugated prior is the Normal-inverse-Wishart, that is basically a distribution over multivariate Normals.

The prior parameters that need to be specified are $\mathbf{\mu}_0$ to define the mean, $\mathbf{\Psi}$ to define the covariance, and two scalar values $\kappa_0$ and $\nu_0$ that I would say define how confident we are on the estimation of the first two parameters respectively.

The updated distribution after observing $n$ samples of a $p$-variate Normal has the form

$$\mathbf{P}(\boldsymbol\mu, \mathbf{\Sigma} | \mathbf{X}) \sim \mathrm{NIW} \left( \frac{\kappa_0\boldsymbol\mu_0+n\mathbf{\bar{x}}}{\kappa_0+n} ,\, \kappa_0+n,\, \nu_0+n ,\, \boldsymbol\Psi + \mathbf{C} + \frac{\kappa_0 n}{\kappa_0+n}(\mathbf{\bar{x}}-\boldsymbol\mu_0)(\mathbf{\bar{x}}-\boldsymbol\mu_0)^T \right)$$

where

$$\mathbf{\bar{x}} = {1 \over n} \sum_{i=0}^{n} \mathbf{x_i}$$

$$\mathbf{C} = \sum_{i=1}^n (\mathbf{x_i} - \mathbf{\bar{x}}) (\mathbf{x_i} - \mathbf{\bar{x}})^T$$

so my desired estimated parameters are

$$E (\boldsymbol\mu | \mathbf{X}) = {{\kappa_0\boldsymbol\mu_0+n\mathbf{\bar{x}}} \over{\kappa_0+n} }$$ $$E (\mathbf{\Sigma} | \mathbf{X}) = \frac{ \boldsymbol\Psi + \mathbf{C} + \frac{\kappa_0 n}{\kappa_0+n}(\mathbf{\bar{x}}-\boldsymbol\mu_0)(\mathbf{\bar{x}}-\boldsymbol\mu_0)^T }{ \nu_0 + n - p - 1}$$

-