Lecture 19 The EM Algorithm

In lecture 5 we found the MLE by solving the likelihood equation, in closed form where this is possible and numerically otherwise. In this lecture we consider models with latent variables or missing data, where the likelihood would be easy to maximise if the missing values were known. The expectation-maximisation (EM) algorithm, introduced in this lecture, makes use of this structure: it alternates between filling in the missing information and re-estimating the parameter values. We illustrate the EM algorithm with the help of a worked example, considering a mixture of two normal distributions. This example also gives the non-identifiable model from lecture 4, to which we return in lecture 20. There we will also prove that each round of the EM algorithm cannot decrease the likelihood, and in lecture 21 we will discuss how to implement the EM algorithm in R.

19.1 Incomplete Data

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

\ell_{c}(\theta)=\sum_{i=1}^{n}\log f_{c}(X_{i},Z_{i};\theta),

and the observed-data log-likelihood is given by

\ell(\theta)=\sum_{i=1}^{n}\log f(X_{i};\theta),

where $f$ is the marginal density of a single observation, given by

f(x;\theta)=\sum_{z}f_{c}(x,z;\theta),

with the integral used instead of the 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.

We will use the following example throughout this lecture and the next.

Example 19.1.

Suppose that a population consists of two groups, with proportion $\pi$ for group 1 and $1-\pi$ for group 2. The measurement has distribution $N(\mu_{1},\sigma^{2})$ in group 1 and $N(\mu_{2},\sigma^{2})$ in group 2, where we assume that we know $\sigma^{2}$ , while $\pi$ , $\mu_{1}$ and $\mu_{2}$ are unknown. 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

f(x;\theta)=\pi\,\varphi_{\sigma}(x-\mu_{1})+(1-\pi)\,\varphi_{\sigma}(x-\mu_{% 2}),

where the parameter is

\theta=(\pi,\mu_{1},\mu_{2})\in(0,1)\times\mathbb{R}\times\mathbb{R}.

This is an example of a two-component normal mixture. The observed-data log-likelihood is

\ell(\theta)=\sum_{i=1}^{n}\log\bigl{(}\pi\,\varphi_{\sigma}(X_{i}-\mu_{1})+(1% -\pi)\,\varphi_{\sigma}(X_{i}-\mu_{2})\bigr{)},

with a sum of two terms inside each logarithm. 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.

19.2 The Algorithm

The EM algorithm starts from a guess $\theta_{0}$ and then alternates between two steps: it replaces $\ell_{c}$ by its conditional expectation given the observed data, computed under the current parameter value, and it maximises the result over $\theta$ .

Definition 19.2.

Let $\theta_{0}\in\Theta$ be given. Then the EM algorithm for finding the parameter values $\theta_{1},\theta_{2},\dots$ consists of the following iterative steps for $k=0,1,2,\dots$ :

(E)

Expectation step: compute

$Q(\theta\mskip 1.0mu|\mskip 1.0mu\theta_{k})=\mathbb{E}_{\theta_{k}}\bigl{(}% \ell_{c}(\theta)\!\mathrel{\big{|}}\!X\bigr{)},$

the conditional expectation of the complete-data log-likelihood given the observed data, where the conditional distribution of $Z$ is computed using the parameter value $\theta_{k}$ .
(M)

Maximisation step: set

$\theta_{k+1}=\mathop{\textnormal{arg\,max}}_{\theta\in\Theta}Q(\theta\mskip 1.% 0mu|\mskip 1.0mu\theta_{k}).$

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

g(z\mskip 1.0mu|\mskip 1.0mux;\theta)=\frac{f_{c}(x,z;\theta)}{f(x;\theta)}.

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

Q(\theta\mskip 1.0mu|\mskip 1.0mu\theta_{k})=\sum_{i=1}^{n}\sum_{z}g(z\mskip 1% .0mu|\mskip 1.0mux_{i};\theta_{k})\,\log f_{c}(x_{i},z;\theta).

For a discrete latent variable the complete-data log-likelihood is linear in the indicators $\mathbf{1}_{\{Z_{i}=z\}}$ , since

\log f_{c}(x_{i},Z_{i};\theta)=\sum_{z}\mathbf{1}_{\{Z_{i}=z\}}\,\log f_{c}(x_% {i},z;\theta).

Comparing this with the formula for $Q(\theta\mskip 1.0mu|\mskip 1.0mu\theta_{k})$ , we see that the E step replaces each indicator $\mathbf{1}_{\{Z_{i}=z\}}$ by its conditional probability $g(z\mskip 1.0mu|\mskip 1.0mux_{i};\theta_{k})$ given the data.

We now apply the algorithm to the two-component normal mixture from example 19.1. We will see that the M step gives the complete-data MLEs from lecture 5, with each group indicator replaced by the conditional probability that the observation belongs to the group.

Example 19.1 (continued).

In this model, the latent variables $Z_{1},\dots,Z_{n}$ are i.i.d. with $\mathbb{P}(Z_{i}=1)=\pi$ and $\mathbb{P}(Z_{i}=2)=1-\pi$ , and given $Z_{i}=j$ we have $X_{i}\sim N(\mu_{j},\sigma^{2})$ . The complete-data density is $f_{c}(x,1;\theta)=\pi\varphi_{\sigma}(x-\mu_{1})$ and $f_{c}(x,2;\theta)=(1-\pi)\varphi_{\sigma}(x-\mu_{2})$ , and thus the complete-data log-likelihood is

\ell_{c}(\theta)=\sum_{i=1}^{n}\Bigl{(}\mathbf{1}_{\{Z_{i}=1\}}\bigl{(}\log\pi% +\log\varphi_{\sigma}(x_{i}-\mu_{1})\bigr{)}+\mathbf{1}_{\{Z_{i}=2\}}\bigl{(}% \log(1-\pi)+\log\varphi_{\sigma}(x_{i}-\mu_{2})\bigr{)}\Bigr{)},

which is linear in the indicators. For the E step, using Bayes’ rule with the current value $\theta_{k}=(\pi_{k},\mu_{1,k},\mu_{2,k})$ , we find

\gamma_{i,k}=\mathbb{P}_{\theta_{k}}(Z_{i}=1\mskip 1.0mu|\mskip 1.0muX_{i}=x_{% i})=\frac{\pi_{k}\,\varphi_{\sigma}(x_{i}-\mu_{1,k})}{\pi_{k}\,\varphi_{\sigma% }(x_{i}-\mu_{1,k})+(1-\pi_{k})\,\varphi_{\sigma}(x_{i}-\mu_{2,k})}.

This is called the responsibility of group 1 for the observation $i$ . By this linearity, taking the conditional expectation replaces $\mathbf{1}_{\{Z_{i}=1\}}$ by $\gamma_{i,k}$ and $\mathbf{1}_{\{Z_{i}=2\}}$ by $1-\gamma_{i,k}$ , and thus

Q(\theta\mskip 1.0mu|\mskip 1.0mu\theta_{k})=\sum_{i=1}^{n}\Bigl{(}\gamma_{i,k% }\bigl{(}\log\pi+\log\varphi_{\sigma}(x_{i}-\mu_{1})\bigr{)}+(1-\gamma_{i,k})% \bigl{(}\log(1-\pi)+\log\varphi_{\sigma}(x_{i}-\mu_{2})\bigr{)}\Bigr{)}.

For the M step we maximise over $\theta=(\pi,\mu_{1},\mu_{2})$ with the $\gamma_{i,k}$ fixed. Since $\log\varphi_{\sigma}(x-\mu)=-(x-\mu)^{2}/(2\sigma^{2})+\text{const}$ , the partial derivatives of $Q$ are

\displaystyle\frac{\partial Q}{\partial\pi}

\displaystyle=\frac{\sum_{i}\gamma_{i,k}}{\pi}-\frac{\sum_{i}(1-\gamma_{i,k})}% {1-\pi},

\displaystyle\frac{\partial Q}{\partial\mu_{1}}

\displaystyle=\frac{1}{\sigma^{2}}\sum_{i=1}^{n}\gamma_{i,k}(x_{i}-\mu_{1}),% \qquad\frac{\partial Q}{\partial\mu_{2}}=\frac{1}{\sigma^{2}}\sum_{i=1}^{n}(1-% \gamma_{i,k})(x_{i}-\mu_{2}),

and setting these to zero gives the updates

equation (19.1) (19.1)

\pi_{k+1}=\frac{1}{n}\sum_{i=1}^{n}\gamma_{i,k},\qquad\mu_{1,k+1}=\frac{\sum_{% i}\gamma_{i,k}x_{i}}{\sum_{i}\gamma_{i,k}},\qquad\mu_{2,k+1}=\frac{\sum_{i}(1-% \gamma_{i,k})x_{i}}{\sum_{i}(1-\gamma_{i,k})}.

The second derivatives are negative and thus these are maxima. They are the complete-data MLEs of examples 5.3 and 5.4 with indicators replaced by responsibilities.

Summary.

•

In an incomplete-data model, the complete data $(X,Z)$ has an easy log-likelihood $\ell_{c}$ , but only $X$ is observed and the observed-data log-likelihood $\ell$ has a sum inside the logarithm.
•

EM iteratively computes the expectation $Q(\theta\mskip 1.0mu|\mskip 1.0mu\theta_{k})=\mathbb{E}_{\theta_{k}}\bigl{(}% \ell_{c}(\theta)\!\mathrel{\big{|}}\!X\bigr{)}$ and then maximises the expectation $Q$ over $\theta$ .
•

For a two-component normal mixture, the E step computes the responsibilities $\gamma_{i,k}=\mathbb{P}_{\theta_{k}}(Z_{i}=1\mskip 1.0mu|\mskip 1.0muX_{i}=x_{% i})$ and the M step computes the weight as the average responsibility and the means as responsibility-weighted means of the data.

Exercise 19.1.

Let $X$ be geometrically distributed with parameter $p\in(0,1)$ , i.e. $f(x;p)=(1-p)^{x-1}p$ for $x\in\mathbb{N}$ , as in exercise 5.1.

1.

Show that $\mathbb{P}_{p}(X\geq m)=(1-p)^{m-1}$ for $m\in\mathbb{N}$ and use this result to show that $\mathbb{P}_{p}(X\geq m+k\mskip 1.0mu|\mskip 1.0muX\geq m)=\mathbb{P}_{p}(X\geq k% +1)$ for all $k\geq 0$ . Comment upon your result in the light of the memorylessness of the exponential distribution from exercise 5.6.
2.

Show that

$\mathbb{E}_{p}(X\mskip 1.0mu|\mskip 1.0muX\geq m)=m-1+\mathbb{E}_{p}(X)=m-1+% \frac{1}{p}\qquad\text{for all $m\in\mathbb{N}$.}$

Exercise 19.2.

Consider the two-component normal mixture from example 19.1, but with known weight $\pi\in(0,1)$ and known variance $\sigma^{2}$ . The only unknown quantities are the means $\theta=(\mu_{1},\mu_{2})$ .

1.

Show that the EM updates for this case are given by the two mean updates from (19.1), where the responsibilities are computed using the known $\pi$ .
2.

Compute the derivative of the observed-data log-likelihood and show that

$\frac{\partial\ell}{\partial\mu_{1}}(\theta)=\frac{1}{\sigma^{2}}\sum_{i=1}^{n% }\gamma_{i}(\theta)\,(x_{i}-\mu_{1}),$

where $\gamma_{i}(\theta)$ is the responsibility at $\theta$ . Deduce that the likelihood equations coincide with the fixed-point equations of the EM iteration.
3.

What are the updates in the limit $\pi\to 1$ ? Comment upon your result in the context of example 5.4.