What is the average bias in MLE of 2-component univariate Gaussian mixture model?

Question

Imagine that you have a standard 2-component univariate Gaussian mixture model:

$$p(x_i∣θ)=(1-λ)N(x_i|μ_1,σ_1^2 )+λN(x_i|μ_2,σ_2^2 )$$ $$θ=\{μ_1,μ_2,σ_1,σ_2,λ\}$$ $$L(θ;x)=∏_{i=1}^N p(x_i |θ)$$

The posterior distribution of parameters can then be obtained using the Bayes rule, and either the mean or the maximum of the posterior can be used as an estimate.

If the parameters are unknown, they can be estimated by standard solutions like MLE using the EM algorithm (which for uniform priors will match the MAP). However, when the number of observations is small and the distributions are not well separated, the estimates will be biased (e.g., https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4257062/). The average bias can of course be estimated through simulations but is there a way to do so analytically for either the maximum or mean of the posterior (with uniform priors)? I am particularly interested in biases for the means $\mu_1$ and $\mu_2$.

Answer 1

There are several issues with mixture model estimates that should be considered.

First, if the variances ($\sigma_{1,2})$ are unknown, then the MLE does not exist, since the likelihood diverges at the points $(\mu_k=x_i,\sigma_k=0)$.

Second, the model has a symmetry $\mu_1\leftrightarrow \mu_2, \sigma_1\leftrightarrow \sigma_2, \lambda \leftrightarrow 1-\lambda$ which makes it non-identifiable. This means, for example, that the posterior probability of the means $P(\mu_1,\mu_2)$ is symmetric $P(\mu_1,\mu_2)=P(\mu_2,\mu_1)$, meaning that $\mathbb E[\mu_1]=\mathbb E[\mu_2]$, namely the posterior means of $\mu_1$ and $\mu_2$ would be the same, and they would lie somewhere in the middle between the two true means, making them not particularly useful estimates by themselves.

Finally as for analytic calculations, since the likelihood is given by

$$\mathcal L = \prod_{i=1}^n \left\{(1-\lambda)f(x_i|\mu_1,\sigma_1) + \lambda f(x_i|\mu_2,\sigma_2) \right\}$$

by expanding the product you can write the posterior distribution itself as a mixture distribution with $2^n$ components, corresponding to all possible partitions of the dataset $\{x_i\}$ to two subsets. Each component would corresponds to independently observing the two subsets from the two gaussian components, and would therefore have a "simple" form:

$$P(\theta|x) = \sum_{\text{partitions}} \text{Beta}(\lambda|n_1+1,n_2+1)P(\mu_1,\sigma_1| \{x_i^{(1)}\}) P(\mu_2,\sigma_2| \{x_i^{(2)}\}) $$

You can calculate various expectations from this distribution analytically as long as you can sum those $2^n$ terms - of course, this becomes impractical very fast.