Consider a mixture of two normal distributions: $ f(x) = p N(x|u_1, S_1) + (1-p) N(x|u_2, S_2) $ where N() is the normal pdf. $p$, $S_2$, and $S_2$ are known. The means are not.
You can get the MLE of $u_1$ and $u_2$ from the EM algorithm.
My question is, for a sample of size n, what are the variances of the EM estimates? My intuition tells me $S_1/(n \times p)$ and $S_2/(n \times (1-p))$ respectively, but I am not sure. Further, do the two estimates have a covariance? What would it be?
I cannot answer with a simple formula to calculate the variances of $\mu_1$ and $\mu_2$ or their covariance. I can only depict the mathematical steps needed to obtain them which leads to the conclusion that there is not an exact analytical solution.
Let's recapitulate for a moment how the variance of an ML estimate $\mu$ for a simple normal distribution with known variance $S$ can be calculated. ML estimation means to select a $\mu$ that maximizes the probability for the data points $\{x_i\}$. Since data points are assumed to be drawn independently, the probability to obtain a sequence of $N$ data points is given by:
$$P(x_1,x_2,\dots)=\prod_{i=1}^{N}\frac{1}{\sqrt{2\pi S}}\exp\left(-\frac{1}{2S}(x_i-\mu)^2\right)=\\ = \left(\frac{1}{\sqrt{2\pi S}}\right)^N \exp\left(-\frac{1}{2S}\sum_{i=1}^N (x_i-\mu)^2\right)$$
In order to select the $\mu$ which maximizes the $P$ one has to find a $\mu$ for which the derivative $\frac{dP}{d\mu}$ vanishes. This can be done analytically and leads to the well known expression $$\mu=\frac{1}{N}\sum_{i=1}^{N}x_i$$ Since the variance $S$ is known, the variance for the sample mean is $$\mbox{Var}(\mu)=\mbox{Var}\left(\frac{1}{N}\sum_{i=1}^{N}x_i\right)=\frac{1}{N}\sum_{i=1}^{N}\mbox{Var}(x_i)=\frac{S}{N}$$
Okay, now for the real objective. The same steps as for the simple normal distribution have to be applied. For the following I abbreviate the density function of the normal distribution by $N$. So, first the probability to obtain a sequence of independently drawn $x_i$ given they follow normal mixture with two components: $$P(x_1,x_2,\dots)=\prod_{i=1}^{N} \left[pN\left(x_i|\mu_1,S_1\right)+(1-p)N\left(x_i|\mu_2,S_2\right)\right]$$ The derivatives of that with respect to $\mu_1$ and $\mu_2$ can be found by applying the chain and the product rule. However, since it is ugly work to do and the result additionally varies depending on the number of data points in a non-trivial way, I cannot state a general analytic expression. Nevertheless, let's assume one calculated the derivatives and obtained two (non-linear) equations to determine $\mu_1$ and $\mu_2$: $$\frac{dP(x_1,x_2,\dots|\mu_1,\mu_2)}{d\mu_1}=0 \hspace{1em}\mbox{and}\hspace{1em} \frac{dP(x_1,x_2,\dots|\mu_1,\mu_2)}{d\mu_2}=0$$ In these equations the unknown variables $\mu_1,\mu_2$ are contained both as factors to be multiplied by the exponential function as well as arguments of the exponential functions. Additionally, the appearing exponential functions have different arguments, so one cannot get rid of them by dividing as it is the case for the simple normal distribution. Hence, these equations cannot be solved analytically. However, if the $\mu_1$ and $\mu_2$ are far apart from each other in terms of $S_1$ and $S_2$ the two stated equations approximately decouple, so that the first equation only contains significant terms with $\mu_1$ and $x_i$'s in vicinity of $\mu_1$ and vice versa for the second equation. Then, the variances of $\mu_1$ and $\mu_2$ will be in good approximation what you already guessed. Furthermore, the correlation between the two variables should be small to zero. However, if there is a significant overlap between the two normal distributions the equations do not decouple.
To conclude, the possibility to state the variance of a quantity explicitly is chained to the possibility to state the quantity itself explicitly. The fact that one has to use an iterative solving scheme such as the EM-algorithm to get $\mu_1,\mu_2$ indicates that one also has to rely on numerical methods for the estimation of their variances.
So, given $S_1,S_2,p$ are known how could one estimate the variance of $\mu_1$ and $\mu_2$? One way is to calculate the two dimensional integral numerically:
$$ \mbox{Var}(\mu_1) = \int_{-\infty}^{+\infty} \int_{-\infty}^{+\infty} (\tilde{\mu}_1 - \mu_1)^2 P(x_1,x_2,\dots|\tilde{\mu}_1,\tilde{\mu}_2) \,d\tilde{\mu}_1\,d\tilde{\mu}_2$$
where $P(x_1,x_2,\dots|\mu_1,\mu_2)$ is the normalized likelihood of the normal mixture:
$$\int_{-\infty}^{+\infty} \int_{-\infty}^{+\infty} P(x_1,x_2,\dots|\tilde{\mu}_1,\tilde{\mu}_2) \,d\tilde{\mu}_1\,d\tilde{\mu}_2 = 1$$
The formula for $\mbox{Var}(\mu_1)$ can be justified within the framework of Bayesian statistics utilizing an improper uniform prior:
$$P(\mu_1,\mu_2|x_1,x_2,\dots) = \frac{P(x_1,x_2,\dots|\mu_1,\mu_2)P(\mu_1,\mu_2)}{P(x_1,x_2,\dots)}$$
The trick is to write down your likelihood exactly. There are many approaches to EM normal mixture estimation. Is the number in each group known? (hypergeometric likelihood) or will you potentially have everyone in one group (binomial likelihood)? Do the 2 normal mixture variates have common variance $\sigma^2$ or are there two separate variances?
Using the law of total variance:
$$\mbox{var}(\bar{x}_1) = \mbox{var}(E(\bar{x}_1|p)) + E(\mbox{var}(\bar{x}_1|p))$$
Where the $\bar{x}_1$ denotes the mean of the first normal mixture variate.
Expressions for the two terms on the RHS depend on the likelihood!
First note that EM is an optimization algorithm. The statistical properties of your estimator are driven by the fact that it is an MLE, and not by the optimization algorithm used to find it.
In this particular case, the MLE has no closed form solution so that you cannot really analyze its properties. Asymptotic theory of MLE does, however, hold for this problem, so that you can use the fact that your estimators are CAN to approximate their properties.
