# Estimate covariance given noisy data and the mean

Basically, I'm inferring the parameters of a Gaussian, $$\mu$$ and $$\Sigma$$ given observed data $$y_i$$ that have uncertainties $$\sigma^2_i$$ associated with them.

# Prior

I intend to use the conjugate prior to the multivariate Gaussian:

• A Gaussian prior on the mean $$\sim \mathcal{N}(\lambda, r)$$

• A Wishart prior on the covariance $$\sim \mathcal{W}(\nu+n, V^{-1})$$

So thats four ($$\lambda, r, \nu, V$$) hyperparameters for the prior.

# Distribution of the mean given noisy data

I can derive the distribution of the mean, $$\mu$$ given the noisy data, $$\{y_i,\sigma^2_i\}$$, which is Gaussian (from the product of gaussian likelihoods):

$$P(\mu|y_i,\sigma_i^2,\Sigma\,\lambda,r) = \mathcal{N}(\mu|\lambda,r)\prod_i \mathcal{N}(y_i|\mu,\Sigma+\sigma^2_i) \propto \mathcal{N}(\mu|m_c, s_c)$$

$$s_c = [(r+\Sigma)^{-1} + \sum_i (\sigma_i^2 + \Sigma)^{-1}]^{-1}$$

$$m_c = s_c [\lambda (r+\Sigma)^{-1} + \sum_i (\sigma_i^2+\Sigma)^{-1}y_i]$$

So I can sample means easily when conditioned on the rest of the parameters.

# Distribution of the covariance given noisy data

Now I need $$P(\Sigma|y_i,\sigma_i^2,\mu\,\nu, V) = \mathcal{W}(\nu+n, V^{-1})\prod_i\mathcal{N}(y_i|\mu,\Sigma+\sigma^2_i)$$

This one is not trival as far as I can tell.

## 1D version

From McLawrence's Question - "estimate-the-variance-of-gaussian-distribution-from-noisy-sample", where they focus on one shared noise variance $$\sigma_e$$ in 1D, the multivariate version will look like the multivariate version of "truncated" inverse gamma distribution (generalises to the truncated inverse Wishart?).

## Version with no data uncertainty

Furthermore, if we say that there are no observational errors ($$\sigma_i=0$$), then we can just use the conjugate expression for Multivariate Normal with known mean found here. This would be an inverse-wishart distribution.

# Conclusion

But, as I say my data are heteroscedastic and so the uncertainties are not the same, making both of the above solutions invalid. We cannot take $$(y - \mu)^T(y-\mu)$$ out of the equation as a constant since there is an extra term $$\sigma_i^2$$ which needs to be included.

Since the distribution for the mean can be readily sampled from, I was wondering if the same can be said for the covariance. Is there a sampling distribution for $$P(\Sigma|y_i,\sigma_i^2,\mu\,\nu, V) = \mathcal{W}(\nu+n, V^{-1})\prod_i\mathcal{N}(y_i|\mu,\Sigma+\sigma^2_i)$$?

For a 1d version this problem I can plot this conditional distribution. It looks like we should be able to sample from it...  In the above plots, the mean is known.

Obviously, this method of plotting cannot be extended to higher dimensions.