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.
 A: $\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.  
A: 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)$$
