Lecture 21 EM in R

Throughout, the parameter is the vector $\theta=(\pi,\mu_{1},\mu_{2})$ and the common standard deviation $\sigma$ is known. The E step of example 19.1 computes the responsibilities $\gamma_{i,k}$ and the M step computes the updates (19.1). Both steps need only vector arithmetic and the normal density, which the command dnorm can be used to compute. To check the iteration, we also compute the observed-data log-likelihood $\ell$, which by proposition 20.2 cannot decrease from one iteration to the next.

We first test the functions on the data of exercise 20.2, the five observations $x=(0,1,5,6,7)$ with $\sigma=1$ and the initial value $\theta_{0}=(1/2,2,4)$. We run one iteration and compare the responsibilities, the iterate $\theta_{1}$ and the two values of the log-likelihood with the hand computation.

Everything agrees with the hand computation ($\mu_{2,1}$ differs in the last digit because the hand computation rounded the intermediate sums).

The two groups in the five observations of exercise 20.2 are far apart, and the iteration converges to $(0.4,0.5,6.0)$, as stated in section 20.2. We now simulate a larger sample in which the groups overlap. For each observation we first choose a group, using the command rbinom to set $Z_{i}=1$ with probability $\pi=0.3$. For the other group R gives the value $0$ instead of $2$, and thus z == 0 marks group 2. We then simulate the observation from the normal distribution of its group, with mean $0$ or $3$ and $\sigma=1$. The command ifelse can be used to choose the mean for each observation.

The box below contains the complete script for the simulation and the EM iteration, with the lines shuffled except for the body of the loop. Before you can run the script, you have to put these lines into a working order. The loop uses the stopping rule of section 20.2: the command break ends it as soon as the log-likelihood increases by less than $10^{-8}$. The vector ll stores the log-likelihood after each iteration, with the initial value at position $1$.

The order follows the rule of lecture 9: each line must come after all lines on which it depends. As in that lecture, several orders satisfy the rule. One of them is the seed, n, z, x, the quartiles q, the start value and the first entry of ll, followed by the loop. The seed and n can also be swapped. The start value puts the weight at $1/2$ and the two means at the lower and upper quartiles, which is a simple way to split the sample into two groups. After the loop, the variable k holds the number of the last iteration.

Before you look at the output, answer two questions. Will the loop need fewer than $10$ iterations, between $10$ and $50$, or more than $50$? Will the log-likelihood increase by the same amount in each iteration?

The loop stops at the fixed point $\hat{\theta}=(0.240,0.220,3.076)$. Its weight is close to the proportion $48/200=0.24$ of group-1 observations which were actually drawn, rather than to $\pi=0.3$, and its two means are within about two tenths of the true values. By proposition 20.3 the fixed point is a root of the likelihood equation.

The command diff can be used to compute the differences of consecutive entries of a vector, and thus d holds the increases of the log-likelihood in the $28$ iterations. The last two lines of the output show the first seven increases and seven increases near the end. They are all positive, as proposition 20.2 guarantees, but they are not equal. Towards the end, each increase is about $0.4$ times the previous one, for example $3.3\cdot 10^{-6}$, then $1.3\cdot 10^{-6}$ and $5.4\cdot 10^{-7}$. This is the geometric convergence mentioned in section 20.2. Since $0.66\cdot 0.4^{20}$ is approximately $7\cdot 10^{-9}$, at this rate about $20$ further iterations are needed to reduce the increase from $0.66$ in the seventh iteration to below the tolerance $10^{-8}$, which explains the total of $28$ iterations. The factor $0.4$ depends on the data, and the take-home task below shows what happens when the groups overlap more.

At the fixed point the responsibility of group 1 is a decreasing function of $x$ which equals $1/2$ where $\pi\varphi_{\sigma}(x-\mu_{1})=(1-\pi)\varphi_{\sigma}(x-\mu_{2})$. Before you look at figure 21.1, predict where the crossing lies: at the midpoint $1.648$ of the two fitted means, or above or below this point?

The crossing lies at $1.244$, well below the midpoint $1.648$, since the weights are unequal: an observation half-way between the means is about three times as likely to come from the larger group as from the smaller one. The filled and open circles give the true group of each observation, which is not known to EM. Far from the crossing point the responsibilities agree with the true groups, but between $0.7$ and $1.6$ observations from both groups occur, and the value of an observation there does not decide its group. The figure is drawn by the script R/S21-mixture.R.

Section 20.2 discusses label switching and the danger of a poor start. In the script R/S21-starts.R the loop is packaged as a function em, which is then run on our sample from three starts: the quartile start, the same start with the two means swapped, and the start $(1/2,3,3)$ with both means equal. Before you look at table 21.1, answer two questions. Where does the start with the swapped means end up? What happens to the responsibilities, and thus to the next update, when $\mu_{1}=\mu_{2}$?

The start with swapped means ends in the mirror image $(0.760,3.076,0.220)$ after the same $28$ iterations and with the same log-likelihood. Both parameter values give the same density, and thus the model is not identifiable in the sense of definition 4.10. The labelling in the start determines which of the two maxima we reach.

When $\mu_{1}=\mu_{2}$, the two normal densities in the responsibility cancel, all $\gamma_{i}$ are equal to $\pi$, and the M step sets both means equal to the sample mean while leaving the weight unchanged. The next step changes nothing, and the loop stops after two iterations at $(1/2,2.3915,2.3915)$. This is effectively a single normal distribution, with a log-likelihood $67$ below the maximum. By proposition 20.3 the point is a root of the likelihood equation. Since the density does not depend on $\pi$ when $\mu_{1}=\mu_{2}$, every point $(\pi,\bar{x},\bar{x})$ is such a root, and the weight $1/2$ is the weight of the start.

The point $(1/2,\bar{x},\bar{x})$ is not a local maximum. The log-likelihood at the nearby start $(1/2,2.39,2.41)$ is already $-421.800$. This is slightly above the value $-421.807$ at the point itself, and figure 21.2 shows that EM started there climbs all the way to the maximum. Neither is it a local minimum, since among the points with $\mu_{1}=\mu_{2}$ the log-likelihood is largest at the sample mean. The point is thus a saddle point, and EM stops there just as it would at a maximum. The figure also shows that near the saddle point the log-likelihood increases only slowly at first: the first iteration from the nearby start gains only $0.07$. With a coarse tolerance such as $0.1$, the loop would have stopped there, far away from the maximum.

The published scripts of this lecture can be used as a starting point for the following task. Choose a weight $\pi$ between $0.2$ and $0.5$ and two means which are two to three units apart. Using your student ID as the seed, simulate a sample of size $n=300$ from the mixture with $\sigma=1$ and run the EM algorithm from the quartile start. The EM algorithm does not provide standard errors. To obtain them, use optim with hessian = TRUE to minimise the negative of loglik, starting from the fixed point, as in lecture 18. The script R/S21-mixture.R contains the necessary command. Report each estimate with its standard error. Then repeat the task with the two means only one unit apart. How many iterations does the loop need now? What happens to the standard errors? If the loop runs through all $1000$ iterations, the stopping rule was never met and the result is not a fixed point. Say so in your answer if this happens.

Typing the commands yourself is an important step towards understanding how the R scripts work. Before the next session, please type in the commands from this session and check that you get the results shown in the text.