EM: the algorithm that alternates guesses

While building a Gaussian Mixture Model from scratch in NumPy for the DLH machine-learning curriculum (the clustering module — same project as my unsupervised clustering notes), I kept tripping over "EM" as a black box. This TIL is the intuition that unlocked it. EM matters beyond the toy: it is the standard way to fit mixture models and latent-variable models — everything from audio segmentation to topic modeling.

You have a cloud of points you believe was generated by k Gaussians — GMM = Gaussian Mixture Model, data modeled as a weighted sum of k bell-shaped distributions (Gaussian = the normal distribution; k = how many clusters you assume) — but you know neither the parameters nor which Gaussian made each point. That's a chicken-and-egg:

• If you knew the model, you'd know the memberships
• If you knew the memberships, you'd know the model

EM breaks the deadlock by alternating guesses:

initialize ──► E-step ──► M-step ──► E-step ──► M-step ──► … until the likelihood stops improving

The three knobs

A GMM models data as a weighted sum of k Gaussians. Each component has three properties:

Identity as a starting covariance is pedagogical: "no correlation between features, equal variance in all directions." EM reshapes it from the data.

The E-step: Bayes in disguise

Bayes' rule answers one question: "given that I saw this data point, which component likely produced it?"

P(k|x) = P(k)·P(x|k) / P(x) — prior × likelihood ÷ total probability. The concrete weather example: 20% of days rainy, 80% sunny; if rainy, 90% chance of clouds, if sunny 20%. You see clouds — what's the chance it's raining?

1. P(clouds) = 0.2·0.9 + 0.8·0.2 = 0.34
2. P(rain|clouds) = 0.2·0.9 / 0.34 ≈ 0.53

53% rain, 47% sun — those must sum to 1. That's exactly the division g = g / column sums. The name "expectation step": g is a soft assignment — every point belongs fractionally to every component instead of being hard-labeled. The log-likelihood must be computed before normalizing (after, every column sum is 1 and log(1) = 0), which is also why the pdf — probability density function, how likely a point is under a Gaussian — floors at 1e-300: a zero pdf would make log(0) = −inf.

The M-step: soft statistics

The M-step does the reverse: fix the responsibilities and ask "what parameters best explain them?" It recomputes each cluster's statistics weighted by the fractions:

1. Soft count: N_k = Σₙ g[k,n] — replaces "number of points" from k-means
2. New prior: π_k = N_k / n — automatically sums to 1
3. New mean: m_k = (Σₙ g[k,n]·xₙ) / N_k — k-means' update, weighted
4. New covariance: S_k = (Σₙ g[k,n]·(xₙ−m_k)(xₙ−m_k)ᵀ) / N_k — the soft scatter matrix

Intuition check: if g were hard (1 for the nearest cluster, 0 otherwise), all four formulas collapse to plain statistics — cluster size ÷ n, cluster mean, cluster covariance. The M-step is just soft statistics. This is the moment the model learns.

The loop is hill-climbing

Each full pass — E then M — is guaranteed to raise the log-likelihood (or keep it flat): fix the parameters → E gives the best possible g; fix g → M gives the best possible parameters. Like a climber who only moves up, holding one coordinate fixed while adjusting the other.

The subtle part: the stop check goes before the M step. ll and g come from the current parameters — break and return immediately and your tuple is self-consistent. Check after the M step and you'd return new parameters with the likelihood of the old ones — a mismatch. The official output shows the climb: -652797 → -94855 → ... → -94439 (early stop at iteration 52). If you ever see it decrease, your E or M step is buggy — a built-in self-test.

Questions that made it click

"How do we apply Bayes' rule here? I'm not following." — Prior P(k) = pi[k]; likelihood P(x|k) = pdf under component k's Gaussian; P(x) = sum over all components of pi[j]·pdf(x|m[j],S[j]) — the law of total probability: the point must have come from some component.

"g[i] = pi[i]·pdf(X, m[i], S[i]) — isn't this prior × posterior?" — Almost — one term is misnamed. pdf is the likelihood, not the posterior. Prior × likelihood = unnormalized posterior; divide by the marginal and then it's the posterior.

"Why the np.sum(np.log(denom))?" — The likelihood of the dataset is a product of n (number of data points) tiny numbers that underflows to 0. Logs turn the product into a sum, and logs of tiny numbers are large-but-finite negatives — nothing underflows.

"What is nk? n×K?" — No — it's Nₖ written without subscripts, the soft count of cluster i. Each point contributes its g[i,n] fraction (0.53 if point n is 53% in cluster i). If assignments were hard, N_k would literally be the number of points in cluster k. That's also why these computations live outside the per-point loop: both steps need the whole g matrix finished — the Bayes denominator is a cross-component quantity, and dividing each point's row while iterating would claim wrong memberships.