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.



Approximate the variance of multiple normal distributions with the same standard deviation
$\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](https://en.wikipedia.org/wiki/Lehmann%E2%80%93Scheff%C3%A9_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. 


