# Approximate the variance of multiple normal distributions with the same standard deviation

Given a number of normal distributions $$N(\mu_1, \sigma^2), N(\mu_2, \sigma^2), ..., N(\mu_n, \sigma^2)$$ with fixed variance $$\sigma^2$$, but not necessary equal means. My question is how to approximate the variance given a number of samples of the normal distributions. Hence given samples $$X^1_1, X^1_2, ..., X^1_{m_1} \sim N(\mu_1, \sigma^2), \\ X^2_1, X^2_2, ..., X^2_{m_2} \sim N(\mu_2, \sigma^2), \\ \vdots \\ X^n_1, X^n_2, ..., X^n_{m_n} \sim N(\mu_n, \sigma^2).$$ Where $$m_1, m_2, ..., m_n > 0$$, but again not necessary equal. What is a good way to approximate the variance $$\sigma^2$$?

With good way to approximate I mean the following. I could take $$m_i$$ such that $$m_i \geq m_1, m_2, ..., m_n$$ and approximate $$\sigma^2$$ with $$\frac{1}{m_i} \sum_{j=1}^{m_i} (X^i_j - E(X^i))^2$$ where $$E(X^i) = \frac{1}{m_i} \sum_k X^i_k$$ is the average over the samples $$X^i_1, ..., X^i_{m_i}$$. Can we do better? A good way to approximate $$\sigma^2$$ would be an estimation method that approximate $$\sigma^2$$ with better precision (on average) than the method described above. I also want to know how one would prove that one estimation method for the variance is better than another method.

For example, my gut feeling is telling me that a weighted average over all variances would approximate better, i.e. $$\frac{1}{m_1 + ... + m_n} \sum_{i=0}^{n} \sum_{j=1}^{m_i} (X^i_j - E(X^i))^2,$$ but I don't know how to prove this. Also I'm worried that one of the variance could be skewed if one of the normal distributions has way less samples than all the others.

$$\newcommand{\si}{\sigma}$$ Let us assume that the $$n$$ samples from the respective distributions $$N(\mu_1, \sigma^2), \dots, N(\mu_n, \sigma^2)$$ are independent. Let $$X_{ij}:=X^i_j$$. Everywhere here $$i=1,\dots,n$$ and $$j=1,\dots,m_i$$. So, all the $$X_{ij}$$'s are independent and $$X_{ij}\sim N(\mu_i, \si^2)$$. So, the joint pdf of $$X:=(X_{ij})$$ is given by $$\begin{multline} f(x)=f_{\mu_1,\dots,\mu_n,\si^2}(x) =(2\pi)^{-m/2}\si^{-m}\exp\Big(-\frac1{2\si^2}\,\sum_{i,j}(x_{ij}-\mu_i)^2\Big) \\ =\exp\Big(-\frac1{2\si^2}\,\sum_{i,j}x_{ij}^2+\sum_i\frac{\mu_i}{\si^2}\,\sum_j x_{ij}\Big) \,c(\mu_1,\dots,\mu_n,\si^2), \end{multline}$$ where $$\begin{equation} m:=\sum_i m_i,\quad x:=(x_{ij}), \end{equation}$$ and $$c(\mu_1,\dots,\mu_n,\si^2)$$ does not depend on $$x$$. So, the $$(n+1)$$-variate statistic $$\begin{equation} S:=\Big(\sum_{i,j}X_{ij}^2,\sum_j X_{1j},\dots,\sum_j X_{nj}\Big) \end{equation}$$ is complete and sufficient for $$(\mu_1,\dots,\mu_n,\si^2)$$.

Let $$\bar X_{i\cdot}:=\frac1{m_i}\,\sum_j X_{ij}$$. Then \begin{align} T&:=\sum_i\Big(\sum_j(X_{ij}-\mu_i)^2-m_i(\bar X_{i\cdot}-\mu_i)^2\Big) \\ &=\sum_i\Big(\sum_j X_{ij}^2-m_i\bar X_{i\cdot}^2\Big) \\ &=\sum_{i,j}X_{ij}^2-\sum_i m_i\bar X_{i\cdot}^2 \end{align} is a function of the complete sufficient statistic $$S$$. Moreover, \begin{align} ET&=\sum_i\Big(\sum_j E(X_{ij}-\mu_i)^2-m_i E(\bar X_{i\cdot}-\mu_i)^2\Big) \\ &=\sum_i(m_i\si^2-m_i \si^2/m_i)=(m-n)\si^2. \end{align} Assume now that $$m>n$$; that is, $$m_i>1$$ for at least one $$i$$. Then the statistic \begin{align} R:=\frac T{m-n}&=\frac1{m-n}\,\sum_i\Big(\sum_j X_{ij}^2-m_i\bar X_{i\cdot}^2\Big) \\ &=\frac1{m-n}\,\sum_{i,j} (X_{ij}-\bar X_{i\cdot})^2 \end{align} is an unbiased estimator of $$\si^2$$, and $$R$$ is also a function of the complete sufficient statistic $$S$$. So, by the Lehmann–Scheffé theorem, $$R$$ is the (essentially unique) uniformly minimum-variance unbiased estimator (UMVUE) of $$\si^2$$.

Notes: (i) I don't think it's a good idea to use the expectation symbol $$E$$ as in your notation $$E(X^i)$$ to denote the arithmetic mean $$\bar X_{i\cdot}$$ of the $$i$$th sample. (ii) With this caveat, the last displayed expression in your post (which is actually the maximum likelihood estimator of $$\si^2$$ here) comes pretty close to the UMVUE $$R$$, except that the factor $$\frac1{m_1 + \dots + m_n}=\frac1m$$ should be replaced by $$\frac1{m-n}$$, to get the unbiasedness; of course, $$\frac1m\sim\frac1{m-n}$$ if $$m$$ is much greater than $$n$$, which should usually be the case.

• Thank you, this is what I had to know. Also thank you for the comments/notes on which symbols to use. – Noud Jan 8 '19 at 18:43

You can solve this precisely computing the basic integrals.

For two normal distributions $$N(\mu_1, \sigma_1)$$ and $$N(\mu_2, \sigma_2)$$ the variance is:

$$\frac{1}{4} \left(\text{\mu 1}^2-2 \text{\mu 1} \text{\mu 2}+\text{\mu 2}^2+2 \text{\sigma 1}^2+2 \text{\sigma 2}^2\right)$$

For three distributions:

$$\frac{1}{9} \left(2 \text{\mu 1}^2-2 \text{\mu 1} \text{\mu 2}-2 \text{\mu 1} \text{\mu 3}+2 \text{\mu 2}^2-2 \text{\mu 2} \text{\mu 3}+2 \text{\mu 3}^2+3 \text{\sigma 1}^2+3 \text{\sigma 2}^2+3 \text{\sigma 3}^2\right)$$

And so forth.

If they have the same standard deviation (or variance):

$$\frac{1}{9} \left(2 \text{\mu 1}^2-2 \text{\mu 1} \text{\mu 2}-2 \text{\mu 1} \text{\mu 3}+2 \text{\mu 2}^2-2 \text{\mu 2} \text{\mu 3}+2 \text{\mu 3}^2+9 \sigma ^2\right)$$