MIT 18.650 — Lecture 12

§ Computing the MLE

In this lecture, we'll turn our attention to actually computing the MLE. For very simple models, computing the MLE as a function of the sample data $\{X_i\}_{i = 1}^n$ is very easy; recall Lecture 9.

More generally, the MLE is the solution to an optimization problem: for what value $\theta$ is $\ell_n(\theta) := \sum_{i = 1}^n \log f_{\theta}(X_i)$ maximized?

There are two natural approaches to solving optimization problems.

Gradient Ascent. Assuming $\ell_n'(\theta)$ is computable, we can start at some guess $\theta_1$ and search for the maximum by picking some constant $\eta \in (0, 1)$ of our choice and updating our guesses $\{\theta_i\}_{i = 1}^{\infty}$ like so:

$$\theta_{i + 1} = \theta_i + \eta \cdot \ell_n'(\theta_i).$$

Newton-Raphson (Root Finding). The idea is to find a root of $\ell_n'(\theta)$ and hope the critical point we find is a global maximum. The Newton-Raphson algorithm updates guesses $\{\theta_i\}_{i = 1}^{\infty}$ for the root like so:

$$\theta_{i + 1} = \theta_i - \frac{\ell'_n(\theta_i)}{\ell''_n(\theta_i)}.$$

Fortunately, when our model is something like $\{\mathcal{N}(\theta, 1) \mid \theta \in \mathbb{R}\}$, the log likelihood $\ell_n(\theta)$ is not very hard to optimize, because its derivative is simple. (This depends on normal distributions nicely having exponential PDFs!)

But there's a special kind of model that is particularly hard to work with: mixture models.

§ Mixture Models

If we measure the heights of a random group of people, a model like $\{\mathcal{N}(\theta, 1) \mid \theta \in \mathbb{R}\}$ won't work perfectly anymore. One reason is that our sample is a mixture of men and women. We might (simplistically) assume the following:

• Any given sample has some probability $p$ of being from a woman, and $(1 - p)$ of being from a man.

• The height of a randomly chosen man has some PDF $f_{\mu_0}$ from the model $\{\mathcal{N}(\mu_0, 1) \mid \mu_0 \in \mathbb{R}\}$.

• The height of a randomly chosen woman has some PDF $f_{\mu_1}$ from the model $\{\mathcal{N}(\mu_1, 1) \mid \mu_1 \in \mathbb{R}\}$.

Our samples thus follow the PDF $f = (1 - p)f_{\mu_0} + pf_{\mu_1}$. This is a mixture model, with parameter $(\mu_0, \mu_1, p)$. More generally…

Definition. In a mixture model, there are two models with parameters $\theta_0$ and $\theta_1$, along with corresponding PDFs $f_{\theta_0}$ and $g_{\theta_1}$. Every sample $i$ yields some data $(Y_i, Z_i)$, where:

$$Z_i \sim \mathrm{Ber}(p) ~~~ \text{ and } ~~~ Y_i \sim \begin{cases} \mathbb{P}_{\theta_0} & \text{ if } Z_i = 0. \\ \mathbb{P}_{\theta_1} & \text{ if } Z_i = 1. \end{cases}$$

Note that we may only receive the data $\{Y_i\}_{i = 1}^n$ and don't even know the value of $p$.

There are two challenges with finding the MLE for mixture models.

Challenge #1. We don't know the value of $p$, so the log likelihood $\ell_n(\mu_0, \mu_1, p)$ has a floating $p$. This means we have a third parameter to worry about now.

Challenge #2. Even if we did assume, say, $p = \frac{1}{2}$, the log likelihood still looks awful.

The difficulty is that the logarithm of a sum does not simplify well at all.

§ The EM Algorithm

It turns out there is a strategy for computing the MLE of a mixture model: the EM Algorithm.

Idea #1. If we know all of the $Z_i$, then estimating $(p, \theta_0, \theta_1)$ is doable.

Reasoning: Of course, if we know all the $Z_i$, we can just estimate $p$ directly.

$$\hat{p} = \frac{1}{n} \sum_{i = 1}^n Z_i.$$

As for $(\theta_0, \theta_1)$, the rough idea is just to split the samples $\{Y_i\}_{i = 1}^n$ into two halves based on their value of $Z_i$. Then just use one half to estimate $\theta_0$, and the other half to estimate $\theta_1$. Formally, this looks like:

$$\hat{\theta}_0 = \argmax_{\theta_0} \sum_{i = 1}^n (1 - Z_i) \log f_{\theta_0}(Y_i) ~~~ \text{ and } ~~~ \hat{\theta}_1 = \argmax_{\theta_1} \sum_{i = 1}^n Z_i \log g_{\theta_1}(Y_i).$$

And these are a lot easier to compute.

Of course, we don't know all of the $Z_i$. But maybe we can at least estimate the $Z_i$…

Idea #2. If we have some prior estimate for $(p, \theta_0, \theta_1)$, then we can estimate $Z_i$ reasonably well.

Reasoning: Well, if we know $(p, \theta_0, \theta_1)$, then $Z_i$ can be estimated via Bayes' Rule.

$$\hat{Z}_i = \frac{g_{\theta_1}(Y_i) \cdot p}{ g_{\theta_1}(Y_i) \cdot p + f_{\theta_0}(Y_i) \cdot (1 - p) }.$$

Here, $\hat{Z}_i$ is not a discrete value in $\{0, 1\}$, but rather a continuous value in $[0, 1]$. Still, it's a good enough estimate that the reasoning in Idea #1 still holds.

We have a problem though: Idea #1 relies on Idea #2, and vice versa! Here's the trick, though: we don't care.

EM Algorithm. Start with some prior estimate $(\theta_0^{(0)}, \theta_1^{(0)}, p^{(0)})$. At time step $t$, make the following updates:

• E Step. In accordance with Idea #2, we estimate $\hat{Z}_i^{(t)}$ using last time step's estimates $(\theta_0^{(t - 1)}, \theta_1^{(t - 1)}, p^{(t - 1)})$ and Bayes' Rule.

  $$\hat{Z}_i^{(t)} = \frac{g_{\theta_1^{(t - 1)}}(Y_i) \cdot p^{(t - 1)}}{ g_{\theta_1^{(t - 1)}}(Y_i) \cdot p^{(t - 1)} + f_{\theta_0^{(t - 1)}}(Y_i) \cdot (1 - p^{(t - 1)}) }.$$

• M Step. In accordance with Idea #1, we compute $(\theta_0^{(t)}, \theta_1^{(t)}, p^{(t)})$ using the estimates $\hat{Z}_i^{(t)}$ computed in the E Step.

  $$p^{(t)} = \frac{1}{n} \sum_{i = 1}^n \hat{Z}_i^{(t)} ~~~ \text{ and } ~~~ \theta_0^{(t)} = \argmax_{\theta_0} \sum_{i = 1}^n (1 - \hat{Z}_i^{(t)}) \log f_{\theta_0}(Y_i) ~~~ \text{ and } ~~~ \theta_1^{(t)} = \argmax_{\theta_1} \sum_{i = 1}^n \hat{Z}_i^{(t)} \log g_{\theta_1}(Y_i).$$

It turns out the following is true.

Proof omitted—it's fairly involved.