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

Question

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.

Answer 1

$\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.

Answer 2

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)$$