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.