ML7.4 · LESSON 4 OF 6

Gaussian mixtures and EM

Perform a clustering update and explain dimensionality reduction, mixtures and anomaly detection.

THE IDEA

When a point could belong to two groups

k-means says every point belongs to exactly one cluster, even a point sitting halfway between two. A Gaussian mixture model (GMM) says the data came from several bell curves mixed together, each with its own mean, spread and weight, and a point belongs to each component with some probability. It does not know which component produced a point, so it computes how responsible each component is.

Example in one dimension: component 1 has weight 0.6 and is centred at 2; component 2 has weight 0.4 and is centred at 6; both have variance 1. For the point x = 3, weight × density gives 0.145 for component 1 and 0.002 for component 2. Divide each by their total, 0.147, and the responsibilities are 0.988 and 0.012. They are non-negative and sum to 1.

Responsibility of component k for point i

r_ik = pi_k * N(x_i | mu_k, Sigma_k)  /  sum over j of  pi_j * N(x_i | mu_j, Sigma_j)

EM: alternate between two guesses

We want the weights, means and covariances that make the observed data most likely, but we cannot solve for them directly because we do not know which component produced which point. Expectation-Maximisation (EM) breaks the knot by alternating:

E-step (expect): with the current parameters, compute every point's responsibilities.
M-step (maximise): with responsibilities as fractional counts, refit each component. Weight = (sum of its responsibilities) / n. Mean = responsibility-weighted average of the points. Covariance = responsibility-weighted average of the outer products of deviations from the new mean.

Compare with k-means: if you replace each row of responsibilities by a hard 1 and 0s, the M-step's mean is exactly the k-means centroid. A GMM is the soft version, and it also learns the shape (covariance) and size (weight) of each group.

Start:

Starting guess. Press Next step: the E-step scores every point under both components. Log-likelihood -212.52.

	weight	mean	var x	cov xy	var y
1	0.5	(2, 2)	1	0	1
2	0.5	(8, 6)	1	0	1

Click a point to see its responsibilities worked out.

Press Next step. Points are coloured by blending the two component colours according to their responsibilities. Ellipses are the 1 and 2 standard deviation contours. Click a point to see its responsibilities worked out.

The log-likelihood only climbs

The score EM improves is the log-likelihood: the sum over points of the log of the mixture density. Each full EM iteration cannot decrease it. In the run above it climbs to -113.05 and the components settle at weights 0.62 and 0.38 (the data was generated with 24 and 16 points, which is 0.6 and 0.4). After the first few steps the gains become tiny, and we stop when they fall below a tolerance.

"Never decreases" is not "reaches the best". EM, like k-means, can stop at a local optimum. Start with both means in the left blob and it converges to a log-likelihood of -117.1, which is 4.05 lower than the -113.05 from a better start. Because the likelihood is not a fraction of anything, you cannot read "−113" as good or bad on its own; only compare fits to the same data.

Predict first. After the E-step a point has responsibilities 0.5 and 0.5. What does that mean?

Chapter context

Unsupervised methods expose structure according to their objectives and assumptions. They do not discover a single objectively correct organization of a dataset; changing distance, dimension or density assumptions can change the result.

TRY IT YOURSELF

Part 1. Do it on the plot

40 points from two overlapping blobs. Use the start buttons, step with Next step, and click points to inspect them.

Start:

Starting guess. Press Next step: the E-step scores every point under both components. Log-likelihood -212.52.

	weight	mean	var x	cov xy	var y
1	0.5	(2, 2)	1	0	1
2	0.5	(8, 6)	1	0	1

Click a point to see its responsibilities worked out.

To do: Press Next step until EM stops (the log-likelihood gain falls below 0.0001)
To do: After the first E-step, click a point whose largest responsibility is below 0.8 in the converged fit (a point two components share)
To do: Converge to the best fit (log-likelihood within 0.01 of the best you can find)
To do: Find an initialisation that converges to a worse fit (log-likelihood more than 1 below the best)

Part 2. Do it in code

responsibilities(x, pi, mu, var) for 1-D data x (length n) and K Gaussians returns an (n, K) array: row i holds pi_k * N(x_i | mu_k, var_k), divided by the row total so it sums to 1. Work in logs and subtract the row maximum before exponentiating, so a point far from every mean does not turn into 0/0.
m_step(x, r) returns (pi, mu, var): Nk is the column sum of r, pi = Nk / n, mu = weighted mean, var = weighted mean of squared distances to the new mu, plus 1e-6 so no variance collapses to 0.
log_likelihood(x, pi, mu, var) returns the sum over points of log of the mixture density, again computed stably.

import numpy as np

def responsibilities(x, pi, mu, var):
    # (n, K) soft assignments; rows sum to 1
    pass

def m_step(x, r):
    # return (pi, mu, var)
    pass

def log_likelihood(x, pi, mu, var):
    # sum_i log sum_k pi_k N(x_i | mu_k, var_k)
    pass

x = np.array([0.0, 5.0, 10.0])
print(responsibilities(x, np.array([0.5, 0.5]), np.array([0.0, 10.0]), np.array([1.0, 1.0])))

Ctrl/⌘ + Enter runs it

Stuck, or something not working? You can come back to it any time.

CHECK YOUR UNDERSTANDING

Must each point belong to exactly one component with probability one?

AYesBNo

Chapter reference: math, worked example & code

This section belongs to the whole chapter, so it is the same on every lesson in it and may cover a different lesson's topic. The code for the experiment you just ran is inside the experiment itself.

Math and notation

J=\sum_i\|x_i-\mu_{c_i}\|^2,\quad \mu_k=\frac1{|C_k|}\sum_{i\in C_k}x_i.

For centered data X, PCA uses leading eigenvectors of $X^TX/(n-1)$ , or an equivalent SVD. Centering and feature scaling are separate choices.

For a Gaussian mixture, EM responsibilities are

r_{ik}=\frac{\pi_k\mathcal N(x_i\mid\mu_k,\Sigma_k)}{\sum_j\pi_j\mathcal N(x_i\mid\mu_j,\Sigma_j)}.

Mixing weights pi are nonnegative and sum to one; each row of responsibilities also sums to one. Weighted sufficient-statistic updates differ from hard k-means assignments.

Worked example

With points [0,2,8,10] and starting centers [0,10], assignments are [0,2] and [8,10]. Updated means are 1 and 9. Squared errors are 1,1,1,1, totaling 4. This example gives one iteration; it does not prove that every initialization finds the best partition.

Reference code

The following example uses the Python standard library and does not call a remote service.

points, centers = [0,2,8,10], [0,10]
groups = [[x for x in points if min(range(2), key=lambda j: (x-centers[j])**2) == i] for i in range(2)]
updated = [sum(group)/len(group) for group in groups]
assert updated == [1,9]
assert sum((x-c)**2 for group,c in zip(groups,updated) for x in group) == 4

Python examples are provided for study and local execution; this page does not execute Python.

Common mistakes

Clusters are objective-dependent. PCA maximizes variance, which may not correspond to useful semantic structure. Anomaly detection is not the same as identifying known fraud classes.