# Estimate the variance of Gaussian distribution from noisy sample

I have measured a large data sample from an underlying Gaussian distribution and want to estimate the variance and its error. However, the measured values are noisy with some Gaussian noise with a standard deviation that is approximately known. How can I estimate the error of the sample variance in this case? First I computed the mean squared error $$\frac{2}{n-1} \sigma^4$$ which does not take into account the noise and seems to be too small.

What I also do not understand is, how one can compute the MSE if one does not know the true distribution, thus don't know the true $\sigma$.

• One aspect that may help, is that you can create a pivotal quantity for this problem. – probabilityislogic Feb 25 '18 at 21:05
• @probabilityislogic Can you elaborate on how you would use that or refer me to some helpful material? – McLawrence Feb 25 '18 at 21:12
• As you've described it the variance of the samples is $\sigma^2+\sigma_\epsilon^2$. So if you know $\sigma_\epsilon^2$, you can subtract that from the variance of the samples. (Estimating the precision of that estimate of $\sigma^2$ is another matter.) If you don't know $\sigma_\epsilon^2$, then you can only estimate the sum of the two variances. One way that allows you to get separate estimates is to take 2 or more estimates/measurements for each sample. Then you can use a mixed model to get the separate estimates. If you know $\sigma_\epsilon^2$, why not consider a bootstrap approach? – JimB Feb 26 '18 at 4:12
• @JimB I need some estimate for the error of my variance. By bootstrapping do you mean something like this: $Var[\hat{\theta}] = \frac{\sum_b^B \hat{\theta}_b - \hat{\bar{\theta}}}{B-1}$? How can I use my knowledge of $\sigma_\epsilon$ in this approach? – McLawrence Feb 26 '18 at 8:02

I would answer this question through the use of bayesian estimation with a "non-informative" prior.

notation and setup

Data model... $(y_i|\mu,\sigma^2)\sim N(\mu,\sigma^2+\sigma_e^2)$ for $i=1,\dots,n$

Prior... $p (\mu,\sigma^2)\propto(\sigma^2+\sigma_e^2)^{-2}$

The prior is not favouring one source. We have $\frac {\sigma^2}{\sigma^2+\sigma_e^2}\sim U (0,1)$ - the variance proportion is uniformly distributed in the prior. Now all this leads to a marginal posterior for $\sigma^2$ of

$$p (\sigma^2|DI)\propto(\sigma^2+\sigma_e^2)^{-(n+1)/2-1}\exp\left(-\frac {(n-1)s^2}{2(\sigma^2+\sigma_e^2)}\right)$$

This is a "truncated" inverse gamma distribution for $(\sigma^2+\sigma_e^2)$ with shape parameter $\frac {n+1}{2}$ and scale parameter $\frac {(n-1)s^2}{2}$ where $s^2=\frac {1}{n-1}\sum_{i=1}^n (y_i-\overline {y})^2$ and $\overline {y}=\frac {1}{n}\sum_{i=1}^n y_i$. The truncation is that $(\sigma^2+\sigma_e^2)>\sigma_e^2$

how to make use of this

First I would check if the truncation matters by calculating the probability that $(\sigma^2+\sigma_e^2)\leq \sigma_e^2$ using the inverse gamma cdf. A quick "eyeball" check can also be that $s^2>>\sigma_e^2$ (i.e. most of the variation in the data does not comes from measurement error).

If the truncation is not important, we can just use the normal inverse gamma. The expected value of the above posterior is $E(\sigma^2|DI)\approx s^2-\sigma_e^2$ and the variance is given as $$var(\sigma^2|DI)=var(\sigma^2+\sigma_e^2|DI)\approx\frac {2}{n-3}\left [E(\sigma^2+\sigma_e^2|DI)\right]^2\approx\frac {2}{n-3}s^4$$. You need to have $n>3$ to apply these formulae.

This is pretty much your expression but with $\sigma^2$ replaced by $s^2$.

If the truncation matters, the you'll need to work in terms of the incomplete gamma function. If you define the function $f (k)=\gamma \left (\frac {n+1}{2}+k, \frac {(n-1)s^2}{2\sigma_e^2}\right)$ then you have....for the expected value... $$E (\sigma^2|DI)=s^2\frac{n-1}{2}\frac {f(-1)}{f (0)}-\sigma_e^2$$ And the variance is given by $$var (\sigma^2|DI)=\left [E(\sigma^2+\sigma_e^2|DI)\right]^2\left [\frac{f (-2)f (0)}{f (-1)^2} - 1\right]$$

• Thanks for your effort! I probably need some time to understand that. My statistics background is not that profound. What exactly means the $DI$ and how could $(\sigma^2 + \sigma_e^2) \leq \sigma_e^2$ be violated if the variances are always positive? – McLawrence Feb 26 '18 at 21:12
• The notation "D" is just shorthand for "the data" and "I" is for "other information" – probabilityislogic Mar 8 '18 at 12:15
• The term $\tau=(\sigma^2+\sigma^2_e)$ will assign a non-zero probability to the event $\tau\leq\sigma^2_e$ if we use the inverse gamma distribution without truncation. However, the formulas are analytic if we do this as an approximation. – probabilityislogic Mar 8 '18 at 12:34

If (using the notation of the answer from @probabilityislogic) $\sigma_e^2$ is known, then the maximum likelihood estimator of $\sigma^2$ is

$$\hat{\sigma}^2={{n-1}\over{n}}s^2-\sigma_e^2$$

(or rather the maximum of this and zero). The estimate of the asymptotic variance of $\hat{\sigma}^2$ is given by

$$\frac{2 (n-1)^2 s^4}{n^3}$$

which for large $n$ is very close to

$$\frac{2s^4}{n}$$

And this is very similar to your result although by $\sigma^2$ you must mean $\sigma^2+\sigma_e^2$. So your formula is consistent with the asymptotic theory. Maybe you could expand on what you mean by "seems to be too small" and why you think that.