1 Week 6: Expectation-Maximisation (EM) algorithm

We motivate the Expectation-Maximisation (EM) algorithm with a classical ${19}^{th}$ century problem!

The first attempt at mixture modelling was by Pearson in 1894 when he was presented by the biologist W.F.R. Weldon with 1000 observations on crabs found in Naples, see [3]. The data are ratios of forehead width to body length, a frequency table is available at
www.maths.uq.edu.au/~gjm/DATA/Crab.dat.

The following graph shows a histogram of the data together with a mixture of two Gaussians fitted.

Unnumbered Figure: Link

The red line gives the mixture density obtained using the EM algorithm and the green lines show the individual Gaussian densities, multiplied by the appropriate weighting. Hence the red line is the sum of the two green lines.
See www.math.mcmaster.ca/peter/mix/demex/excrabs.html for more details.

We shall consider the two-group case although the ideas readily extend to more than two groups, and from mixtures of Gaussian to other distributions.

If a random variable $X$ is drawn from $N({\mu }_1,{\sigma }_1^2)$ with probability ${\alpha }_1$ and from $N({\mu }_2,{\sigma }_2^2)$ with probability ${\alpha }_2$, ${\alpha }_1+{\alpha }_2=1$, then the pdf of $X$ is

$$f(x)={\alpha }_1\varphi (x|{\mu }_1,{\sigma }_1^2)+{\alpha }_2\varphi (x|{\mu }_2,{\sigma }_2^2)$$

where $\varphi (x;\mu ,{\sigma }^2)$ is the pdf of $N(\mu ,{\sigma }^2)$. This pdf is called a Gaussian mixture with two components.

Proof. Let $Y=j$ if $X$ come from $N({\mu }_j,{\sigma }_j^2)$, $j=1,2$. By definition, the cdf of $X$ is

	$F(x)$	$=$	$P(X\le x)$
		$=$	$P(X\le x|Y=1)P(Y=1)+P(X\le x|Y=2)P(Y=2)$
		$=$	${\displaystyle {\int }_{-\mathrm{\infty }}^x}\varphi (x;{\mu }_1,{\sigma }_1^2)dx\times {\alpha }_1+{\displaystyle {\int }_{-\mathrm{\infty }}^x}\varphi (x;{\mu }_2,{\sigma }_2^2)dx\times {\alpha }_2$

Thus the pdf

$$f(x)=F^{\prime }(x)=\varphi (x|{\mu }_1,{\sigma }_1^2)\times {\alpha }_1+\varphi (x|{\mu }_2,{\sigma }_2^2)\times {\alpha }_2$$

Given a random sample $x_1,\mathrm{\dots },x_n$ from this distribution, the log-likelihood function is

	$\mathrm{\ell }(\theta ;𝐱)$	$=$	${\displaystyle \sum _{i=1}^n}\log \left\{{\alpha }_1\varphi (x_i|{\mu }_1,{\sigma }_1^2)+{\alpha }_2\varphi (x_i|{\mu }_2,{\sigma }_2^2)\right\}$
		$=$	${\displaystyle \sum _{i=1}^n}\log \left\{{\displaystyle \frac{{\alpha }_1}{\sqrt{2\pi }{\sigma }_1}}{e}^{-\frac{1}{2}{(\frac{x_i-{\mu }_1}{{\sigma }_1})}^2}+{\displaystyle \frac{{\alpha }_2}{\sqrt{2\pi }{\sigma }_2}}{e}^{-\frac{1}{2}{(\frac{x_i-{\mu }_2}{{\sigma }_2})}^2}\right\}$

which cannot be maximised easily. Here $𝐱=(x_1,\mathrm{\dots },x_n)$ and $\theta =({\alpha }_1,{\alpha }_2,{\mu }_1,{\mu }_2,{\sigma }_1^2,{\sigma }_2^2)$.

However, with additional information $y_1,\mathrm{\dots },y_n$ on which distribution each $x_i$ comes from, the likelihood function becomes

	$L(\theta ;𝐱,𝐲)$	$=$	${\displaystyle \prod _{i=1}^n}\mathbb{P}(X_i=x_i,Y_i=y_i|\theta )$
		$=$	${\displaystyle \prod _{i=1}^n}\mathbb{P}(X_i=x_i|Y_i=y_i,\theta )\mathbb{P}(Y_i=y_i|\theta )$
		$=$	${\displaystyle \prod _{i=1}^n}\left\{{\displaystyle \frac{1}{\sqrt{2\pi }{\sigma }_{y_i}}}\exp \left(-{\displaystyle \frac{1}{2}}{\left({\displaystyle \frac{x_i-{\mu }_{y_i}}{{\sigma }_{y_i}}}\right)}^2\right)\times {\alpha }_{y_i}\right\}.$

Thus the log-likelihood function becomes

	$\mathrm{\ell }(\theta ;𝐱,𝐲)$	$=$	${\displaystyle \sum _{i=1}^n}\log \left\{{\displaystyle \frac{1}{\sqrt{2\pi }{\sigma }_{y_i}}}{e}^{-\frac{1}{2}{\left(\frac{x_i-{\mu }_{y_i}}{{\sigma }_{y_i}}\right)}^2}\times {\alpha }_{y_i}\right\}$
		$=$	$-{\displaystyle \frac{n}{2}}\log (2\pi )-{\displaystyle \frac{1}{2}}{\displaystyle \sum _{i=1}^n}\log ({\sigma }_{y_i}^2)-{\displaystyle \frac{1}{2}}{\displaystyle \sum _{i=1}^n}{\left({\displaystyle \frac{x_i-{\mu }_{y_i}}{{\sigma }_{y_i}}}\right)}^2+{\displaystyle \sum _{i=1}^n}\log {\alpha }_{y_i}$
		$=$	$-{\displaystyle \frac{n}{2}}\log (2\pi )-{\displaystyle \frac{1}{2}}\log ({\sigma }_1^2){\displaystyle \sum _{i=1}^n}{1}_{\{y_i=1\}}-{\displaystyle \frac{1}{2}}\log ({\sigma }_2^2){\displaystyle \sum _{i=1}^n}{1}_{\{y_i=2\}}$
			$-{\displaystyle \frac{1}{2}}{\displaystyle \sum _{i=1}^n}{\left({\displaystyle \frac{x_i-{\mu }_1}{{\sigma }_1}}\right)}^2{1}_{\{y_i=1\}}-{\displaystyle \frac{1}{2}}{\displaystyle \sum _{i=1}^n}{\left({\displaystyle \frac{x_i-{\mu }_2}{{\sigma }_2}}\right)}^2{1}_{\{y_i=2\}}$
			$+\log {\alpha }_1{\displaystyle \sum _{i=1}^n}{1}_{\{y_i=1\}}+\log {\alpha }_2{\displaystyle \sum _{i=1}^n}{1}_{\{y_i=2\}}$

which is simpler and can be maximised as follows with ‘obvious’ solutions.

Setting the partial derivatives with respect to ${\mu }_1$ and ${\sigma }_1^2$ to zero

$$\frac{\partial }{\partial {\mu }_1}\mathrm{\ell }(\theta ;𝐱,𝐲)=-\frac{1}{2}\sum _{i=1}^{n}2\left(\frac{x_i-{\mu }_1}{{\sigma }_1}\right)\frac{-1}{{\sigma }_1}{1}_{\{y_i=1\}}=0,$$

$$\frac{\partial }{\partial {\sigma }_1^2}\mathrm{\ell }(\theta ;𝐱,𝐲)=-\frac{1}{2}\frac{1}{{\sigma }_1^2}\sum _{i=1}^n{1}_{\{y_i=1\}}+\frac{1}{2}\frac{1}{{\sigma }_1^4}\sum _{i=1}^n{(x_i-{\mu }_1)}^2{1}_{\{y_i=1\}}=0,$$

gives

$${\widehat{\mu }}_1=\frac{{\sum }_{i=1}^n{x}_i{1}_{\{y_i=1\}}}{{\sum }_{i=1}^n{1}_{\{y_i=1\}}},{\widehat{\sigma }}_1^2=\frac{{\sum }_{i=1}^n{(x_i-{\widehat{\mu }}_1)}^2{1}_{\{y_i=1\}}}{{\sum }_{i=1}^n{1}_{\{y_i=1\}}}$$

which are respectively the sample mean and sample variance of the sample from $N({\mu }_1,{\sigma }_1^2)$. Similarly,

$${\widehat{\mu }}_2=\frac{{\sum }_{i=1}^n{x}_i{1}_{\{y_i=2\}}}{{\sum }_{i=1}^n{1}_{\{y_i=2\}}},{\widehat{\sigma }}_2^2=\frac{{\sum }_{i=1}^n{(x_i-{\widehat{\mu }}_2)}^2{1}_{\{y_i=2\}}}{{\sum }_{i=1}^n{1}_{\{y_i=2\}}}$$

The rest of the maximisation uses the Lagrange multiplier method

$$\log {\alpha }_1\sum _{i=1}^n{1}_{\{y_i=1\}}+\log {\alpha }_2\sum _{i=1}^n{1}_{\{y_i=2\}}-\lambda ({\alpha }_1+{\alpha }_2-1)$$

with constraint ${\alpha }_1+{\alpha }_2=1$.

$$\frac{1}{{\alpha }_1}\sum _{i=1}^n{1}_{\{y_i=1\}}-\lambda =0,\frac{1}{{\alpha }_2}\sum _{i=1}^n{1}_{\{y_i=2\}}-\lambda =0$$

$${\widehat{\alpha }}_1=\frac{1}{n}\sum _{i=1}^n{1}_{\{y_i=1\}},{\widehat{\alpha }}_2=\frac{1}{n}\sum _{i=1}^n{1}_{\{y_i=2\}}$$

which are the sample proportions from each population.

However, the data $𝐲$ are not available in practice. Therefore how can we make use of the above maximum likelihood estimator for $\theta$ given $𝐱$ and $𝐲$ when the only data available are $𝐱$. This is where the EM-algorithm comes into play. We iterate between the $M$-step (MLE) above computing $\widehat{\theta }$ based on the full data ($𝐱$ and $𝐲$) and an $E$-step which estimates the unobserved $𝐲$ given the data $𝐱$ and current parameter estimates $\theta$.

Before formally introducing the EM algorithm, we outline the procedure for the mixture model.

To initialise, we need to choose initial estimates for the parameters of the model. The key thing is to choose ${\mu }_1\ne {\mu }_2$, and for the parameters to be reasonable. For example, the data ranges from 0.57 to 0.7, therefore I choose ${\mu }_1=0.6$ and ${\mu }_2=0.65$ (plausible values in the range) with ${\sigma }_1^2={\sigma }_2^2=0.0004$ (close to the variance of the data $0.000364$). Initially I have no preference for either of the two Gaussian distributions making up the mixture, so take ${\alpha }_1={\alpha }_2=0.5$.

E-Step
The augmented data that we want are $𝐲$, the mixture component (Gaussian distribution) to which each of the observations $𝐱$ belongs. Given $\theta =({\mu }_1,{\mu }_2,{\sigma }_1,{\sigma }_2,{\alpha }_1,{\alpha }_2)$ and $𝐱$, we can compute the probability (expectation) of the observations coming from each of the two distributions.

Since the observations are assumed to be independent,

$$\mathbb{P}(𝐘=𝐲|\theta ,𝐱)=\prod _{i=1}^n\mathbb{P}(Y_i=y_i|\theta ,x_i).$$

By Bayes’ theorem, for $i=1,2,\mathrm{\dots },n$,

	$𝔼[1_{\{Y_i=1\}}|\theta ,x_i]=\mathbb{P}(Y_i=1|\theta ,x_i)$	$=$	${\displaystyle \frac{\mathbb{P}(X_i=x_i,Y_i=1|\theta )}{\mathbb{P}(X_i=x_i|\theta )}}$
		$=$	${\displaystyle \frac{\mathbb{P}(X_i=x_i,Y_i=1|\theta )}{\mathbb{P}(X_i=x_i,Y_i=1|\theta )+\mathbb{P}(X_i=x_i,Y_i=2|\theta )}}$
		$=$	${\displaystyle \frac{{\alpha }_1/(\sqrt{2\pi }{\sigma }_1)\exp \{-{(x_i-{\mu }_1)}^2/(2{\sigma }_1^2)\}}{{\alpha }_1/(\sqrt{2\pi }{\sigma }_1)\exp \{-{(x_i-{\mu }_1)}^2/(2{\sigma }_1^2)\}+{\alpha }_2/(\sqrt{2\pi }{\sigma }_2)\exp \{-{(x_i-{\mu }_2)}^2/(2{\sigma }_2^2)\}}},$

with $𝔼[1_{\{Y_i=2\}}|\theta ,x_i]=1-𝔼[1_{\{Y_i=1\}}|\theta ,x_i]$.

Thus in the E-Step we compute the expectation of $𝔼[1_{\{Y_i=1\}}|\theta ,x_i]$, which we shall call $p_i$ (with $𝔼[1_{\{Y_i=2\}}|\theta ,x_i]=1-p_i$).

M-Step
We can now plug in the expected values for $1_{\{y_i=1\}}$ and $1_{\{y_i=2\}}$, obtained in the E-step into the full data MLEs derived above to update our estimates of $\theta$. Specifically,

The algorithm then iterates between the two steps (E and M). Updating the probability of each observation belonging to each of the two Gaussian distributions and then updating the estimates (via MLEs) of the parameters. The process stops when we have convergence. That is, two successive estimates of the parameters agree to some pre-defined precision.

For the Naples bay crab data I required successive values of ${\alpha }_1$ to differ by less than ${10}^{-4}$. This meant that with my chosen starting values 467 iterations were required (taking 0.22 seconds in R!), resulting in ${\mu }_1=0.6319$, ${\mu }_2=0.6550$, ${\sigma }_1^2=0.000329$, ${\sigma }_2^2=0.000155$, ${\alpha }_1=0.444$ and ${\alpha }_2=0.556$.

The EM algorithm is an iterative algorithm, introduced in [2], and is designed to compute the maximum likelihood estimate when there is missing data. It is iterative in that it consists of a series of steps called iterations, where the parameter value gets repeatedly updated, until a convergence criterion is met. The algorithm converges to a local maximum of the likelihood function (ie. not necessarily where the global maximum is). Thus if the likelihood function is unimodal the EM algorithm will converge to the maximum likelihood estimate (MLE).

The EM algorithm is primarily used for maximising the likelihood when the task becomes easier given more information associated with existing data. This situation is called an incomplete data problem because we do not have this extra information. A prime example is mixture distributions such as that given above. In the Naples crab example, if we knew from which of the two Gaussian distributions each observation comes parameter estimation is straightforward.

The EM algorithm gets its name from the two steps involved in each iteration. The E-step (expectation) and the M-step (maximisation). The procedure alternates between these two steps (so it goes E, M, E, M,…) until convergence is achieved (successive values are sufficiently close).

Let $x$ represent our observed data and $\theta$ the parameter(s) of interest. (In the Naples crab example, $x$ represents the crab ratios and $\theta$ represents the parameters of the two Gaussian distributions and the mixture proportions.) By choosing a suitable parametric model we can write down the likelihood function $L(\theta ;x)$ and our aim is then to find the value of $\theta$ that maximises $L(\theta ;x)$ or equivalently $\mathrm{\ell }(\theta ;x)=\log (L(\theta ;x))$, i.e. find the MLE (maximum likelihood estimate) of $\theta$. Suppose that it is possible to think of some additional data $y$, such that the ‘complete data’ log-likelihood $\mathrm{\ell }(\theta ;x,y)$ is of a simpler form than the ‘incomplete data’ log-likelihood $\mathrm{\ell }(\theta ;x)$ and thus easier to maximise. (In the Naples crab example, $y$ represents from which of the two Gaussian distributions each observation comes.) The log-likelihood of $x$ and $y$ is preferred but $y$ is not available, so we estimate the log-likelihood by taking expectation with respect to the available data $x$ in such a way that the estimated log-likelihood function remains easier to estimate. However, the expectation requires the value of $\theta$, thus it can only be done iteratively, using the current value of $\theta$. These are the key ideas behind the EM algorithm. Starting with an initial value for $\theta$, we

1. 1.

   Estimate the log-likelihood $\mathrm{\ell }(\theta ;x,y)$ using the current value of $\theta$ by taking expectation conditioning on $x$.

2. 2.

   Maximise the estimated log-likelihood with respect to $\theta$ to get an updated value for $\theta$.

3. 3.

   Continue steps 1 and 2 until convergence is achieved.

Note that the above is given in terms of the log-likelihood because it is often easier to work with than the likelihood, of the complete data $x$ and $y$.

More specifically, the procedure is as follows. Let ${\theta }^{*}$ denote the current estimate of $\theta$, we define the function

$Q(\theta ,{\theta }^{*})={\displaystyle \int \mathrm{\ell }(\theta ;x,y)f(y|x;{\theta }^{*})dy}$

that is, in the expectation step,

$$Q(\theta ,{\theta }^{*})=E_{{\theta }^{*}}[\mathrm{\ell }(\theta ;x,Y)|X=x].$$

The expectation is with respect to the conditional distribution of $Y$ given the observed data $x$ using the current estimate ${\theta }^{*}$ of $\theta$, ie. $f(y|x;{\theta }^{*})$. Note that $Q(\theta ,{\theta }^{*})$ is the expectation of the complete data log-likelihood function taking the unavailable data $Y$ as random. The E- and M- steps can then formally be specified as:

1. 1.

   The E-step Calculation of (key elements of ) $Q(\theta ,{\theta }^{*})$ as a function of $\theta$;
   (Determination of $Q(\theta ,{\theta }^{*})$ by taking expectation might sound daunting but it is often fairly straightforward. Only the key elements of $Q$ are calculated, so that it can be maximised in the M-step.)

2. 2.

   The M-step Maximisation of $Q(\theta ,{\theta }^{*})$ with respect to $\theta$ to get ${\theta }^{\ast \ast }$, which becomes ${\theta }^{*}$ in the next iteration.
   (Maximisation of $Q(\theta ,{\theta }^{*})$ should be much easier than maximising $\mathrm{\ell }(\theta ;x)$, often with an explicit solution.)

The EM algorithm can be interpreted as ’data augmentation’ when the complete data log-likelihood is linear in the missing data $y$, which is estimated in the E-step and substituted in the M-step. That is, rather than the somewhat off-putting $Q(\theta ,{\theta }^{*})$, the E-step reduces down to compute $E[y|x,{\theta }^{\ast }]$ and this is then substituted in for $y$ in the M-step. In the Naples crab example, this is not quite the case but we simply computed the probability that each observation belongs to each of the two Gaussian distributions and substitute this via $𝔼[1_{\{y_i=1\}}|x_i,\theta ]=\mathbb{P}(y_i=1|x_i,\theta )$ into the MLE. In other words, we can formally show that the straightforward procedure we employed for the Naples crab data is in fact an implementation of the EM algorithm.

The EM algorithm is therefore a deterministic algorithm that alternates between the expectation and maximisation steps. Even though it works with the (estimated) complete data log-likelihood, it actually increases the incomplete data log-likelihood $\mathrm{\ell }(\theta ;x)$ in each iteration. Let us prove this step by step.

First, we give an overview of the notation. Below we assume the data are continuous for simplicity of exposition but the same arguments (with minor modifications) apply to discrete data or a mixture of continuous and discrete data. Note that the Naples crab data is a mixture of continuous and discrete data, the crab ratios $𝐱$ are continuous, whilst the allocation variable $𝐲$ are discrete. The full data likelihood is $f(x,y|\theta )$, the probability density function of the observed and augmented data given parameters $\theta$. The observed data likelihood is $f(x|\theta )=L(\theta ;x)={\int }_{y}f(x,y|\theta )𝑑y$. Finally, $f(y|x,\theta )$ is the conditional distribution of $y$ given $x$ and $\theta$.

Because we maximise $Q(\theta ,{\theta }^{*})$ as a function of its first argument to get ${\theta }^{**}$,

	$Q({\theta }^{**},{\theta }^{*})-Q({\theta }^{*},{\theta }^{*})$	$\ge$	$0$
	${\displaystyle \int \log (f(x,y|{\theta }^{**}))f(y|x,{\theta }^{*})dy}-{\displaystyle \int \log (f(x,y|{\theta }^{*}))f(y|x,{\theta }^{*})dy}$	$\ge$	$0$
	${\displaystyle \int \log \left(\frac{f(x,y|{\theta }^{**})}{f(x,y|{\theta }^{*})}\right)f(y|x,{\theta }^{*})dy}$	$\ge$	$0.$		(1.1)

Since $f(x,y|\theta )=f(y|x,\theta )f(x|\theta )$, we can write

$$\int \log \left(\frac{f(y|x,{\theta }^{**})}{f(y|x,{\theta }^{*})}\right)f(y|x,{\theta }^{*})dy+\int \log \left(\frac{f(x|{\theta }^{**})}{f(x|{\theta }^{*})}\right)f(y|x,{\theta }^{*})dy\ge 0.$$

(1.2)

Using the inequality $\log (x)\le x-1$, we can deduce that the first integral cannot be positive.

	${\displaystyle \int \log \left(\frac{f(y|x,{\theta }^{**})}{f(y|x,{\theta }^{*})}\right)f(y|x,{\theta }^{*})dy}$	$\le$	${\displaystyle \int \left(\frac{f(y|x,{\theta }^{**})}{f(y|x,{\theta }^{*})}-1\right)f(y|x,{\theta }^{*})dy}$		(1.3)
		$=$	${\displaystyle \int f(y|x,{\theta }^{**})dy}-{\displaystyle \int f(y|x,{\theta }^{*})dy}$
		$=$	$0.$

Therefore the second integral must be non-negative.

	${\displaystyle \int \log \left(\frac{f(x|{\theta }^{**})}{f(x|{\theta }^{*})}\right)f(y|x,{\theta }^{*})dy}$	$\ge$	$0$
	$\log \left({\displaystyle \frac{f(x|{\theta }^{**})}{f(x|{\theta }^{*})}}\right){\displaystyle \int f(y|x,{\theta }^{*})dy}$	$\ge$	$0$
	$\log (f(x|{\theta }^{**}))-\log (f(x|{\theta }^{*}))$	$\ge$	$0.$		(1.4)

Since $f(x|\theta )=L(\theta ;x)$, we have that $\mathrm{\ell }({\theta }^{**};x)\ge \mathrm{\ell }({\theta }^{*};x)$, the log-likelihood is increasing (or at the very least non-decreasing) from one iteration to the next.

We have shown that at each iteration of the EM algorithm the likelihood is increasing. Furthermore, it can be shown that if these iterations converge, then they converge to a stationary point of the likelihood function. Although this is necessarily a local maximum, unless the likelihood function is unimodal this is not necessarily a global maximum, ie. the MLE.

Advantages

1. 1.

   It often has meaningful solutions and interpretation as data augmentation.

2. 2.

   Convergence is easy to ‘detect’ since we wait until successive values of $\theta$ fall within a certain level of variation.

Disadvantages

1. 1.

   It can be slow to converge.

2. 2.

   Need to be able to compute the E-step – restricts applications. (This can be circumvented to some extent by using a Monte Carlo EM algorithm.)

This example actually originates from the paper by [2] which introduced the EM algorithm.

The data consist of the genetic linkage of 197 animals and the data are divided into four genetic categories, labeled 1 through to 4, see [4] and [2]. The probabilities that an animal belongs to each of the four categories are $\frac{1}{2}+\frac{\theta }{4},\frac{1-\theta }{4},\frac{1-\theta }{4},\frac{\theta }{4}$, respectively. Let $𝐱=(x_1,x_2,x_3,x_4)=(125,18,20,34)$ denote the total number of animals in each category. For example, the probability of belonging to category 1 is $\frac{1}{2}+\frac{\theta }{4}$ and there are 125 observed animals in category 1. It is possible to maximise this multinomial likelihood directly and hence, obtain the MLE for $\theta$ without recourse to the EM algorithm. However, it is far simpler to apply the EM algorithm. To show that the EM algorithm works, we do both for this example.

MLE - observed data

The likelihood for $\theta$ given $𝐱$ is

	$L(\theta ;𝐱)$	$=$	${\displaystyle \frac{(x_1+x_2+x_3+x_4)!}{x_1!x_2!x_3!x_4!}}{\left({\displaystyle \frac{2+\theta }{4}}\right)}^{x_1}{\left({\displaystyle \frac{1-\theta }{4}}\right)}^{x_2+x_3}{\left({\displaystyle \frac{\theta }{4}}\right)}^{x_4}$		(1.5)
		$\propto$	${(2+\theta )}^{x_1}{(1-\theta )}^{x_2+x_3}{\theta }^{x_4}.$

Therefore the log-likelihood satisfies

$\mathrm{\ell }(\theta ;𝐱)=K+x_1\log (2+\theta )+(x_2+x_3)\log (1-\theta )+x_4\log (\theta ),$

(1.6)

where $K$ is a constant. Therefore

${\displaystyle \frac{d}{d\theta }}\mathrm{\ell }(\theta ;𝐱)={\displaystyle \frac{x_1}{2+\theta }}-{\displaystyle \frac{x_2+x_3}{1-\theta }}+{\displaystyle \frac{x_4}{\theta }}.$

(1.7)

By setting $\frac{d}{d\theta }\mathrm{\ell }(\theta ;𝐱)=0$, we obtain the following quadratic (in $\theta$) equation from (1.7),

	$x_1(1-\theta )\theta -(x_2+x_3)(2+\theta )\theta +x_4(2+\theta )(1-\theta )$	$=$	$0$
	$-\{x_1+x_2+x_3+x_4\}{\theta }^2+\{x_1-2(x_2+x_3)-x_4\}\theta +2x_4$	$=$	$0$
	$-197{\theta }^2+15\theta +68$	$=$	$0$		(1.8)

This yields

$$\theta =\frac{-15\pm \sqrt{{15}^2-4\times 68\times -197}}{2\times -197}=-0.5507\text{ or }0.6268.$$

Since $\theta$ needs to lie between 0 and 1, we get $\theta =0.6268$.

EM algorithm - full data

What extra information is going to make computation of the MLE simpler?

Suppose that we can choose $y$ so that we have an augmented data likelihood of the form

$L(\theta ;𝐱,y)=C{\theta }^A{(1-\theta )}^B.$

(1.9)

Then

$\mathrm{\ell }(\theta ;𝐱,y)=\log C+A\log \theta +B\log (1-\theta )$

with

${\displaystyle \frac{d}{d\theta }}\mathrm{\ell }(\theta ;𝐱,y)={\displaystyle \frac{A}{\theta }}-{\displaystyle \frac{B}{1-\theta }}.$

Setting $\frac{d}{d\theta }\mathrm{\ell }(\theta ;𝐱)=0$ give the linear equation

$$A(1-\theta )-B\theta =0.$$

It is then trivial to show that

$$\theta =\frac{A}{A+B}.$$

Note that likelihoods of the form (1.9) are common in statistics, for example, binomial (geometric, negative binomial) data.

What is stopping the genetic linkage data from having a likelihood of the form (1.9)? The first cell with probability $1/2+\theta /4$.

What if we break the probability $1/2+\theta /4$ into $1/2$ and $\theta /4$, a component independent of $\theta$ and a component proportional to $\theta$? Below we detail how this can be done.

Suppose that the observed cell $x_1$ could be divided into two subcategories $A$ and $B$. Suppose that $y$ $(x_1-y)$ is the number of animals is subcategory $A$ $(B)$ with cell probability $\theta /4$ $(1/2)$. This would give an augmented data set $(𝐱,y)$.

The likelihood for $\theta$ given $(𝐱,y)$, is

	$L(\theta ;𝐱,y)$	$=$	${\displaystyle \frac{(x_1+x_2+x_3+x_4)!}{y!(x_1-y)!x_2!x_3!x_4!}}{\left({\displaystyle \frac{1}{2}}\right)}^{x_1-y}{\theta }^y{(1-\theta )}^{x_2+x_3}{\theta }^{x_4}$		(1.10)
		$\propto$	${(1-\theta )}^{x_2+x_3}{\theta }^{y+x_4}$

Hence, there exists a constant $C$ such that

	$\mathrm{\ell }(\theta ;𝐱,y)$	$=$	$\log (L(\theta ;𝐱,y))$		(1.11)
		$=$	$\log C+(x_2+x_3)\log (1-\theta )+(y+x_4)\log (\theta ).$

(This yields $\widehat{\theta }=(y+x_4)/(y+x_2+x_3+x_4)$ using the observations after (1.9).)

E - step
From (1.11), we have that

	$Q(\theta ,{\theta }^{*})$	$=$	$E[\mathrm{\ell }(\theta ;𝐱,y)]$
		$=$	$E[\log C|{\theta }^{*},𝐱]+E[(y+x_4)\log (\theta )+(x_2+x_3)\log (1-\theta )|{\theta }^{*},𝐱]$
		$=$	$K+(E[y|{\theta }^{*},𝐱]+x_4)\log (\theta )+(x_2+x_3)\log (1-\theta ),$

since $(x_2,x_3,x_4)$ are known, fixed constants and $K$ does not depend on $\theta$, so plays no role in the M-step (can be ignored). Therefore the only quantity that needs to be computed in the E-step is $E[y|{\theta }^{*},𝐱]$.

What is the distribution of $y|{\theta }^{*},𝐱$?
There are 125 animals in category 1 which independently can be assigned to either category $A$ or category $B$. The conditional probability of animal belonging to category $A$ given that they belong to category 1 is simply

$$\frac{P(\text{Category }A\text{ and Category 1})}{P(\text{Category 1})}=\frac{P(\text{Category }A)}{P(\text{Category 1})}=\frac{{\theta }^{*}/4}{{\theta }^{*}/4+1/2}=\frac{{\theta }^{*}}{{\theta }^{*}+2}.$$

Thus $y|{\theta }^{*},𝐱\sim \mathrm{Bin}(125,{\theta }^{*}/(2+{\theta }^{*}))$, and hence, $E[y|{\theta }^{*},𝐱]=125{\theta }^{*}/(2+{\theta }^{*})$. This gives

$Q(\theta ,{\theta }^{*})=C+\left\{{\displaystyle \frac{125{\theta }^{*}}{2+{\theta }^{*}}}+34\right\}\log (\theta )+38\log (1-\theta ).$

(1.12)

M - step
The M-step is straightforward and simply involves finding $\frac{d}{d\theta }Q(\theta ,{\theta }^{*})=0$ (which we have effectively already done above). From (1.12), we have that

$$\frac{d}{d\theta }Q(\theta ,{\theta }^{*})=\left\{\frac{125{\theta }^{*}}{2+{\theta }^{*}}+34\right\}\frac{1}{\theta }-\frac{38}{1-\theta }.$$

Hence ${\theta }^{\ast \ast }$ solves

$$\left\{\frac{125{\theta }^{*}}{2+{\theta }^{*}}+34\right\}\frac{1}{{\theta }^{\ast \ast }}-\frac{38}{1-{\theta }^{\ast \ast }}=0,$$

which following the steps above yields

$${\theta }^{\ast \ast }=\frac{125{\theta }^{*}/(2+{\theta }^{*})+34}{125{\theta }^{*}/(2+{\theta }^{*})+34+38}.$$

For this example, we have that the maximum likelihood estimate (MLE) is $\theta =0.6268$ to 4 decimal places. Starting with ${\theta }^0=0.5$ as the initial value for $\theta$, successive iterations of the EM algorithm yields:

The EM algorithm is particularly useful where we have censored data. In many experiments and medical trials, we have censored data. In a medical study, we could be interested in the time until death following major surgery. However, we may only be able to follow up patients for 2 years after surgery. That is, for a patient alive 2 years after surgery we don’t know when the death of the patient occurs, only that death occurs more than 2 years after surgery.

Suppose that $X_1,X_2,\mathrm{\dots },X_n$ are independent and identically distributed according to $f(x|\theta )$. However, assume that there is right censoring at the value $a$. That is, for all observations greater than $a$, the only information we have is that the observation exceeds $a$. The common notation is to denote censored values superscripted with an asterisk. Therefore, we have an incomplete sample (suitably rearranged)

$$𝐱=(x_1,\mathrm{\dots },x_m,a^{*},\mathrm{\dots },a^{*}).$$

In this case, the likelihood based on all we know about $x_1,\mathrm{\dots },x_n$ is

$$L(\theta ;𝐱)=\left\{\prod _{i=1}^{m}f(x_i|\theta )\right\}{\{1-F(a|\theta )\}}^{n-m},$$

where $F$ is the cdf (cumulative distribution function) corresponding to $f(x)$. Typically it is difficult to maximise $L(\theta ;𝐱)$ because of the presence of the non-standard term ${\{1-F(a|\theta )\}}^{n-m}$. Let $𝐲=(x_{m+1},x_{m+2},\mathrm{\dots },x_n)$ represent the true but unobserved values. Then the pair $(𝐱,𝐲)$ provides the complete data $(x_1,\mathrm{\dots },x_n)$ and

$$L(\theta ;𝐱,𝐲)=\prod _{i=1}^{n}f(x_i|\theta )$$

which is usually straightforward to maximise.

Suppose that $X\sim Gam(2,\delta )$, where $\delta$ is unknown. Suppose that we have the following 15 observations with the data censored at $a=2.5$:

	$1.226,{2.500}^{\ast },1.229,0.576,1.925$
	${2.500}^{\ast },1.437,1.217,1.836,{2.500}^{\ast }$
	${2.500}^{\ast },1.643,2.225,{2.500}^{\ast },2.067$

Note that the cdf of $X$ is given by

$$F(x)=1-e^{-\delta x}(1+\delta x)x>0,$$

whilst the pdf of $X$ is given by $f(x)={\delta }^{2}x{e}^{-\delta x}$, $x>0$. We reorder the data, $𝐱$, into ascending order, then

$$L(\delta ;𝐱)=\left\{\prod _{i=1}^{10}{\delta }^2{x}_i{e}^{-\delta x_i}\right\}{\left(e^{-2.5\delta }(1+2.5\delta )\right)}^5.$$

This is difficult to maximize.

Let $y_i$ denote the actual value of the $i^{th}$ observation. Then for $i=1,2,\mathrm{\dots },10$, $y_i=x_i$ and for $i=11,\mathrm{\dots },15$, $y_i>x_i$.

Then

	$\mathrm{\ell }(\delta ;𝐱,𝐲)$	$=$	${\displaystyle \sum _{i=1}^{15}}\log \left({\delta }^2{y}_i{e}^{-\delta y_i}\right)$
		$=$	${\displaystyle \sum _{i=1}^{15}}\{2\log (\delta )+\log (y_i)-\delta y_i\}$
		$=$	$30\log (\delta )-\delta {\displaystyle \sum _{i=1}^{15}}{y}_i+C,$

where $C$ does not depend on $\delta$.

Hence,

$$\frac{\mathrm{d}}{\mathrm{d}\delta }\mathrm{\ell }(\delta ;𝐱,𝐲)=\frac{30}{\delta }-\sum _{i=1}^{15}{y}_i$$

which when set equal to 0 solves to give $\widehat{\delta }=2/\overline{y}$, where $\overline{y}$ is the sample mean.

Now conditioning on what we know about $x_1,\mathrm{\dots },x_n$,

$$\mathrm{E}[\mathrm{\ell }(\delta ;𝐱,𝐘)|𝐱]=30\log (\delta )-\delta \left(x_1+\mathrm{\cdots }+x_{10}+5x^{\ast }\right)+C$$

where $x^{\ast }=𝔼[X|X>a]$ using the conditional density of $X$ given $X>a$. This leads to

$$Q(\delta ,{\delta }^{*})=30\log (\delta )-\delta \left(x_1+\mathrm{\cdots }+x_{10}+5x^{\ast }\right)+\text{const}$$

which is maximised at

$${\delta }^{**}=\frac{30}{x_1+\mathrm{\cdots }+x_{10}+5x^{\ast }}.$$

where

$$x^{*}=\frac{2+2{\delta }^{*}a+{\delta }^{*2}{a}^2}{{\delta }^{*}(1+{\delta }^{*}a)}.$$

The EM algorithm consists of calculation of $x^{*}$ in the E-step and ${\delta }^{**}$ in the M-step, which are repeated until convergence.

For the given data, $x_1+\mathrm{\cdots }+x_{10}=15.381$ and the EM algorithm converges to $0.8387$.

Note that

	$x^{\ast }$	$=$	$𝔼[X|X>a]$
		$=$	${\displaystyle {\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}}x{f}_X(x|x>a)𝑑x.$

Now

$$f_X(x|x>a)=\frac{1_{\{x>a\}}{f}_X(x)}{\mathbb{P}(X>a)},$$

where $\mathbb{P}(X>a)=1-F(a)$. Therefore

	$x^{\ast }$	$=$	${\displaystyle {\int }_a^{\mathrm{\infty }}}x{\displaystyle \frac{f_X(x)}{1-F(a)}}𝑑x$
		$=$	${\displaystyle \frac{1}{1-F(a)}}{\displaystyle {\int }_a^{\mathrm{\infty }}}xf(x)𝑑x$
		$=$	${\displaystyle \frac{1}{(1+\delta a)\exp (-\delta a)}}{\displaystyle {\int }_a^{\mathrm{\infty }}}x{\delta }^{2}x\exp (-\delta x)𝑑x$
		$=$	${\displaystyle \frac{1}{(1+\delta a)\exp (-\delta a)}}\times {\displaystyle \frac{1}{\delta }}(2+2\delta a+{\delta }^2{a}^2)\exp (-\delta a)$
		$=$	${\displaystyle \frac{2+2\delta a+{\delta }^2{a}^2}{\delta (1+\delta a)}}.$