Lecture 19 The EM Algorithm

Assume that we are given data consisting of i.i.d. pairs $(X_{i},Z_{i})$, $i=1,\dots,n$, with joint density or joint probability weights $f_{c}(x,z;\theta)$, but we have observed only $X=(X_{1},\dots,X_{n})$. The pair $(X,Z)$ with $Z=(Z_{1},\dots,Z_{n})$ is called the complete data, the values $X$ are called the observed data and the unobserved values $Z_{i}$ are called latent variables. The complete-data log-likelihood is given by

and the observed-data log-likelihood is given by the marginal density of a single observation,

where the integral is used instead of a sum when $Z_{i}$ is continuous. The MLE for $\theta$ is the maximiser of $\ell$, but the algebra from lecture 5 is not applicable, since the sum inside the logarithm prevents any simplifications. The function $\ell_{c}$, being a sum of standard log-likelihoods, would be easy to maximise, but it cannot even be evaluated as it depends on $Z$. The EM algorithm is based on this observation.

In the example we consider in the lecture, the population consists of two groups, with proportion $\pi$ for group 1 and $1-\pi$ for group 2. The measurement for group 1 is $N(\mu_{1},\sigma^{2})$ and for group 2 is $N(\mu_{2},\sigma^{2})$, where we assume that we know $\sigma^{2}$. We are given observations of the measurement $X_{i}$, but we do not know which group $Z_{i}\in\{1,2\}$ the observation belongs to. If we write $\varphi_{\sigma}$ for the density of $N(0,\sigma^{2})$, then the observations are i.i.d. with density

This is an example of a two-component normal mixture. If we knew the groups, the MLEs for this example would be easy to find: we would get the proportion of group 1 and the two group means by using examples 5.3 and 5.4, but without the groups we cannot find the MLEs in closed form.

If we knew the parameter, we could compute the conditional distribution of the latent variables given the data; if we knew the latent variables, we could maximise $\ell_{c}$. Neither of these quantities is known, so we start with a guess $\theta_{0}$ and alternate between filling in the missing values, replacing $\ell_{c}$ by its conditional expectation given the observed data, and maximising over $\theta$.

The second argument $\theta_{k}$ of $Q$ enters only through the conditional distribution of $Z$ given $X$. By Bayes’ rule, proposition A.20 in appendix A, this conditional distribution has density or probability weights

Since the pairs $(X_{i},Z_{i})$ are independent, once the data are observed $Q(\theta\mskip 1.0mu|\mskip 1.0mu\theta_{k})$ is a sum over the observations, given by

In cases where $\ell_{c}$ is linear in the indicator functions of the latent variables, e.g. for the mixture distribution considered in this section, the E-step is to replace the indicator functions by the corresponding conditional probabilities given the data.

The algorithm only ever maximises $Q$, and it is not obvious that this does anything for the observed-data log-likelihood $\ell$. The central result of this lecture is that it does. We assume throughout that the set of values $z$ with $f_{c}(x,z;\theta)>0$ does not depend on $\theta$, the analogue for the latent variables of condition (R2) in definition 16.1. We also assume that $f(x;\theta)>0$ for every data point $x$ and every $\theta\in\Theta$, and that $\mathbb{E}_{\theta^{\prime}}\bigl{(}\log f_{c}(x,Z;\theta)\!\mathrel{\big{|}}\!X=x\bigr{)}$ is finite for all $\theta,\theta^{\prime}\in\Theta$ and every data point $x$, so that $Q(\theta\mskip 1.0mu|\mskip 1.0mu\theta_{k})$ and $H(\theta\mskip 1.0mu|\mskip 1.0mu\theta_{k})$ in the decomposition (19.1) below are well defined. Rearranging the formula for $k(z\mskip 1.0mu|\mskip 1.0mux;\theta)$ and taking logarithms, we have

for every such $z$. The left-hand side does not involve $z$ and thus taking conditional expectations given $X=x$ under $\theta_{k}$, with the random variable $Z$ in place of $z$, does not change this. Summing over the sample gives the decomposition

valid for all $\theta,\theta_{k}\in\Theta$.

The statement of the theorem does not imply that EM finds the MLE: the sequence $\ell(\theta_{k})$ converges whenever $\ell$ is bounded above, but the iteration may stop at a local maximum or, in principle, at a saddle point. In general we can only say that the points where the EM algorithm converges are stationary points of $\ell$.

For a vector parameter the same argument applies, but now using gradients instead of derivatives. Conversely, if the function $\theta\mapsto Q(\theta\mskip 1.0mu|\mskip 1.0mu\theta^{\prime})$ is strictly concave for all $\theta^{\prime}$, then all roots of the likelihood equation are fixed points: The derivative of $\theta\mapsto Q(\theta\mskip 1.0mu|\mskip 1.0mu\theta^{*})$ is zero at a root $\theta^{*}$ and the only maxima of a strictly concave function are the points where the derivative is zero; exercises 19.2 and 19.4 show that this argument applies for two mixtures.

We can apply the algorithm for finding the parameters of a mixture from section 19.1. The E and M steps for this mixture model are as follows.

In the example we assumed that the variance was known and common to both groups. If, on the other hand, we had assumed separate, unknown variances for the two groups, then the estimate for one component could be centred on an observation and the variance of this component could be shrunk to zero. In this case, the likelihood would become infinite at this point and the MLE would not exist.

Three issues arise whenever EM is used, and the mixture shows all three. The first is the promised failure of identifiability: the parameter values $(\pi,\mu_{1},\mu_{2})$ and $(1-\pi,\mu_{2},\mu_{1})$ give the same density $f(x;\theta)$, and thus the model is not identifiable in the sense of definition 4.10, the difficulty foreshadowed after example 4.11. This is called label switching: the likelihood has two global maxima which are mirror images of each other, and EM started from $(1/2,4,2)$ instead of $(1/2,2,4)$ in example 19.4 converges to $(0.6,6.0,0.5)$ instead of $(0.4,0.5,6.0)$. Either answer is fine on its own, but averaging different runs is not, and the remedy of lecture 4 is to restrict the parameter space by a convention such as $\mu_{1}<\mu_{2}$. The second issue is local maxima: theorem 19.2 guarantees ascent and nothing more, and a poor starting value can lead the iteration to a local maximum of $\ell$ without warning. The standard precaution is to run EM from several starting values and to keep the run with the largest observed-data log-likelihood, which we can always evaluate. The third issue is when to stop: since $\ell(\theta_{k})$ is non-decreasing, we stop when its increase falls below a small tolerance; convergence is geometric. Finally, EM produces an estimate but no standard error; that is obtained afterwards from the observed information at the fixed point, as in section 17.3. Lecture 21 shows all three issues in R.