Gaussian Mixtures and Anomaly Detection: Clustering With Shape

Lectures 11a and 11b. A direct continuation of the K-means post: what if the groups aren't round?

Where K-means limps again

The professor generates a dataset on purpose: two stretched, rotated groups (multiplying the points by a rotation matrix), plus a much smaller third group set apart from the other two. K-means, even with a good initialization (centers picked by hand, close to the right spot), struggles: since it only ever sees "distance to a center," it tends to cut the elongated groups into rounder pieces than they really are, because its boundary between two groups is always a straight line perpendicular to the line joining the centers, never an ellipse.

Gaussian mixtures: every group becomes an ellipse, not a point

from sklearn.mixture import GaussianMixture
gm = GaussianMixture(n_components=3, n_init=10, random_state=42)
gm.fit(X)

The core difference from K-means: instead of storing just one center per group, a Gaussian mixture stores an entire normal distribution per group, with its own mean and covariance matrix (which captures the shape, how elongated the group is and in which direction). The full model is a weighted sum of $K$ Gaussians:

$p(\mathbf{x}) = \sum_{k=1}^{K} \pi_k \, \mathcal{N}(\mathbf{x} \mid \boldsymbol{\mu}_k, \boldsymbol{\Sigma}_k)$

where $\pi_k$ is the weight (the fraction of data belonging to that group). Bishop derives fitting this model via an algorithm called EM (expectation-maximization), which alternates two steps until it converges:

• E step (expectation): for every point, compute each component's responsibility $\gamma(z_{nk})$, the probability (via Bayes) that the point came from that specific group, given where the Gaussians currently sit. Unlike K-means, which assigns each point to one group only, EM assigns a fractional responsibility across all groups (a point on the border between two groups might get 60% responsibility from one and 40% from the other).
• M step (maximization): recompute each Gaussian's mean, covariance, and weight, using those responsibilities as weights. A point with 0.9 responsibility toward group 1 counts almost fully toward group 1's mean and covariance, while one split 0.5/0.5 counts half toward each.

This is literally the "soft" version of K-means: swap "each point belongs to exactly one group" for "each point belongs a little to every group," and swap "a group is just a mean" for "a group is a mean plus a shape."

Output: weights found, [0.40, 0.21, 0.39] (matches the true proportion of the three generated groups). Converged in just 4 iterations.

Interactive: the ellipses settling into place

My own reconstruction of the EM algorithm (on the same 380 points, 3 groups, two of them stretched). Click "EM step" and watch the ellipses (each the 1-standard-deviation outline of that Gaussian) rotate and stretch until they fit the data's real shape:

Notice the weights shown below the chart drift toward 0.40, 0.21, 0.39 as you click "EM step" repeatedly, and the ellipses go from generic circles (the initialization) to stretched shapes that trace the groups' real direction.

Anomaly detection for free

One advantage of having a density model (not just a grouping): you can ask "how likely is this point, given the model?" and flag the least likely ones as anomalies.

densities = gm.score_samples(X)
density_threshold = np.percentile(densities, 2)
anomalies = X[densities < density_threshold]

Points with density below the 2nd percentile (the least likely 2%) become anomaly candidates, and they capture exactly the small isolated group the professor placed on purpose far from the other two.

How many groups to use? BIC, AIC, and a smarter way

gm.bic(X), gm.aic(X)

Output: BIC = 8189.73, AIC = 8102.51.

Both are information criteria: they measure how well the model explains the data, with a penalty for complexity (more groups = more parameters = bigger penalty), to avoid picking "more groups is always better" just because more groups always fits better. Running for $K=1$ through $9$ and plotting both against $K$, the shape of the curve points at the number of groups that balances fit against simplicity.

There's an even more direct way: BayesianGaussianMixture is handed a generous number of components (10, in this case) and prunes the unnecessary ones itself, zeroing out their weight:

from sklearn.mixture import BayesianGaussianMixture
bgm = BayesianGaussianMixture(n_components=10, n_init=10, random_state=42)
bgm.fit(X)
print(np.round(bgm.weights_, 2))

Output: [0.4, 0.21, 0.39, 0, 0, 0, 0, 0, 0, 0]. Seven of the ten components zeroed out on their own, leaving the three real ones, no need for me to scan $K$ by hand.

Bishop warns of a technical problem worth knowing: if a Gaussian "collapses" right on top of a single data point, its variance can go to zero and the model's likelihood goes to infinity, a singularity, not a good fit. Real implementations (like scikit-learn's) guard against this in practice with numerical safeguards, but it's a reminder that "finding the maximum likelihood" isn't always as well-behaved a problem as it sounds.

Switching topics: other ways to spot an anomaly

aula11b sets up a more direct anomaly-detection scenario: 980 "normal" points (3 well-behaved groups) plus 20 points scattered randomly across the space (the real anomalies), and compares three different detectors, each with the default contamination (0.1, meaning "assume 10% of the data is anomalous") and then tuned via Optuna:

The default contamination=0.1 tells every detector to flag 10% of the data as anomalous (98 points), but only 20 of the 1000 points (2%) are actually anomalies. Forcing the model to find 5 times more anomalies than actually exist guarantees a pile of false positives, hence the low F1 across all three before tuning. Once Optuna searches for the right contamination (close to 0.02, the true value) and each model's other hyperparameters, all three improve substantially, and LocalOutlierFactor (which decides "anomalous" by comparing a point's local density against its neighbors', the same density-based reasoning DBSCAN uses) comes out ahead.

The practical lesson, matching the rest of the lecture: contamination isn't a cosmetic detail, it's the most important piece of the tuning, because it tells the model how many anomalies to look for. Without knowing (or estimating well) that fraction ahead of time, any of these three detectors misses the mark.

Wrapping up

Practical application

I use the same density idea from the GMM anomaly-detection section, but compare four different percentile cutoffs on the same 3-group dataset, to see how much the cutoff choice changes how many points turn into "anomalies."

for percentil in [1, 2, 5, 10]:
    threshold = np.percentile(densities, percentil)
    n_anomalias = (densities < threshold).sum()
    print(percentil, n_anomalias)

Output (1250 points total, the same 750+250 generated in the lecture): with a 1% cutoff, 13 points flagged. With 2% (the value used in the notebook), 25 points. With 5%, 63 points. With 10%, 125 points.

The number of "anomalies" found scales almost linearly with the chosen percentile, because that's literally how a percentile cutoff works: it always finds exactly that fraction of the data, whether a real anomaly is sitting there or not. It's the same point the lecture's detector comparison already made a different way: deciding "how much" to look for is a choice that changes the result as much as the algorithm itself.