Maximum likelihood estimation and bayesian inference of variance given multiple datasets

Question

I'm currently working on a problem were I have multiple normal distributed data sets $X_1, \dotsc,X_n$ with each data set having it's own mean $\bar x_i $ but all have the same variance $\sigma$. The sample size differs between the sets and can also reach 1 for some $X_i$. Now I want to estimate $\sigma^2 $ ($\bar x_i $'s are unknown). What would be the best way of doing this? $\bar\sigma^2 =1/(N-1)\sum_{i=1}^n\sum_{x\in X_i} (x-\bar x_i)^2$, where $N$ is total number of samples. This would underestimate the variance at least in the case their is a 1 sample data set $X_i$ as then $(x-\bar x_i)^2=0$, which would give no information about the variance.

Should the normalization look some think like $1/(N-n)$?
Furthermore, if I want to expand this to a Bayesian inference approach, how would the loglikelihoodfunction look like. The default I think would look like
$\sum_{i=1}^n\sum_{x\in X_i} \frac{(x-\bar x_i)^2}{\sigma^2}-N*log(\sigma\sqrt{2\pi})$,
which would presumably also underestimate $\sigma$. My idea would be replace $N$ with $\hat N$, which is the total number of samples in data sets with more than 1 sample. Is this sufficient?

Edit:
the true means of $X_i$ are all unknown but in the case of MLE are estimated by $\bar x_i=1/n\sum_{x\in X_i}x$. (should still be the MLE) In the case of Bayesian inference I also want to handle them as a unknown variable, i.e. $\bar x_i$ should also be derived via the inference process.

Answer 1

Under your assumptions (multiple samples from populations with the same variance $\sigma^2$) the optimal estimator of $\sigma^2$ is $$ \frac{1}{n - k} \sum_{j=1}^{k} \sum_{i=1}^{n_j} (x_{ij} - \bar{x}_{j})^2 = \frac{1}{n - k} \left[ (n_1-1)s_1^2 + (n_2-1)s_2^{2} + \ldots + (n_{k}-1) s_k^2 \right] $$ where for sample $j = 1,2, \ldots, k$ of size $n_j$, each $s_j^2 = \frac{1}{n_j -1} \sum_{i=1}^{n_j} (x_{ij} - \bar{x}_{j})^2$ is a sample estimate of $\sigma^2$ and thus the whole sum above is a weighted average across your samples (which, by the way, explains $\frac{1}{n-k} = \frac{1}{\sum_{j} (n_j - 1)} $).

You might recognize this expression as the mean sum of squares due to the error from the ANOVA F-ratio, where it provides exactly what you're asking for: a pooled estimate of a common variance parameter.