Gaussian mixture models

What this chapter covers

k-means gave every point one hard label and drew every cluster as a round ball. That worked on four tidy blobs, and it falls apart the moment the groups are stretched, tilted, or overlapping — which is most of the time. This chapter keeps the good idea from k-means, that a cluster has a center, and fixes the two things that hurt it. A cluster gets to be a tilted, stretched ellipse instead of a circle, and a point gets to belong partly to more than one of them instead of being forced to pick.

That model is the Gaussian mixture, and we fit it with expectation-maximization — the same soft-assign-then-refit loop from the EM chapter, now in two dimensions with a full covariance matrix per component so each cluster can carry its own shape and angle. We build the whole thing by hand in NumPy: the multivariate Gaussian, the E-step that hands out soft responsibilities, the M-step that re-fits every mean, covariance, and weight, and the log-likelihood that the loop climbs. Then we let scikit-learn's GaussianMixture fit the same model and check that our number lands on theirs.

The centerpiece is an animation of the ellipses molding themselves to the data — three round blobs that peel apart and tilt into place, points shading from one color into another where the clusters overlap, a log-likelihood panel rising underneath. And the payoff is a side-by-side you can't unsee: the same data clustered by hard spherical k-means and by soft elliptical GMM, one of them slicing straight across groups the other one wraps cleanly.

A bit of history

The pieces are old. In 1894 Karl Pearson, working on a set of crab measurements that came out lopsided instead of bell-shaped, decided the sample was really two species mixed together and set out to recover both bells from the blend. He did it by the method of moments — matching five sample moments to the parameters of a two-component Gaussian mixture, which left him solving a ninth-degree polynomial by hand. It's one of the first times anyone fit a mixture model to real data, and he did it without a computer, without maximum likelihood, and without any of the machinery that would later make it routine.

The machinery arrived in 1977, when Arthur Dempster, Nan Laird, and Donald Rubin named the expectation-maximization algorithm and showed that a huge family of "the data I have is an incomplete version of data I could fit easily" problems all reduced to the same loop. A Gaussian mixture is the textbook case: if you knew which component made each point, fitting the Gaussians would be one weighted average; if you knew the Gaussians, labeling the points would be one Bayes rule. You know neither, so EM alternates between them. Pearson's crabs finally had a general method, and the mixture of Gaussians became the standard first example of EM in every course since — including the last one.

The intuition

Start with what k-means can't do. It measures everything by plain distance to a center, so its idea of a cluster is a ball: equally wide in every direction, the same size everywhere, boundaries that are straight lines halfway between centers. Hand it a cluster shaped like a long diagonal streak and it has no way to say so. It plants a center in the middle and calls the far ends outliers, or worse, gives them to a neighbor whose center happens to sit closer.

A Gaussian mixture replaces the ball with a full Gaussian. A Gaussian in two dimensions has a mean, which is the center, and a covariance matrix, which is the shape: how wide, how tall, and crucially how tilted. That covariance is the whole upgrade. A single variance would give you a circle back; a diagonal covariance gives you an axis-aligned ellipse; a full covariance, with the off-diagonal term free, gives you an ellipse rotated to any angle. Each cluster wraps itself in the ellipse that actually fits its points.

The second upgrade is softness. Instead of assigning a point to one cluster, the model asks: given the current ellipses, what's the probability this point came from each one? A point deep inside one ellipse gets a near-certain answer. A point sitting in the overlap between two gets a split — 60% this one, 40% that one — and the model keeps that ambiguity instead of pretending it away. Those probabilities are called responsibilities, and they're both how the model fits and what it hands back.

Here's the data we'll work on. Three blobs, but not round ones: they've been run through a shear that stretches and tilts them into long diagonal ellipses that lean the same way and overlap at the seams. No labels, no colors — just 300 points on a plane.

Your eye still finds three groups, but notice they're diagonal and they touch. That's the data k-means was built to fail on and the mixture was built to handle.

The math

A mixture says the data comes from $K$ Gaussians blended together. Each component $k$ has a mixing weight $\pi_k$ (its share of the data, with the weights summing to one), a mean vector $\mu_k$, and a covariance matrix $\Sigma_k$. The probability the model assigns to a point $x$ is the weighted sum of what each component thinks of it:

The shape term is the multivariate Gaussian density. For a $D$-dimensional point it is:

The quantity $(x-\mu)^{\top}\Sigma^{-1}(x-\mu)$ is the squared Mahalanobis distance — ordinary squared distance warped by the inverse covariance, so that "far" is measured in units of the ellipse's own spread and direction. When $\Sigma$ is the identity this collapses to Euclidean distance and the level sets are circles, which is exactly the k-means view. A full $\Sigma$ is what tilts and stretches those circles into ellipses.

Fitting means choosing all the $\pi_k, \mu_k, \Sigma_k$ to make the data as likely as possible. EM does it in two steps. The E-step computes the responsibility of component $k$ for point $i$ — the posterior probability it came from that component, which is just its weighted density normalized across components:

Every row $\gamma_{i\cdot}$ sums to one; that's the soft assignment. The M-step then re-fits each component as a responsibility-weighted version of the ordinary sample estimates. Let $N_k = \sum_i \gamma_{ik}$ be the soft count of points in component $k$. The new weight is its share:

The new mean is the weighted average of the points:

And the new covariance is the weighted scatter around that mean — the equation that lets each ellipse take its own shape and tilt:

Each point contributes to each component in proportion to how much it belongs there. The number the loop is climbing is the total log-likelihood — the log of the mixture probability, summed over points:

EM's guarantee is that $\ell$ never decreases from one iteration to the next. When it stops climbing, we stop.

What it's good at, what it isn't

The mixture wins wherever the round-ball assumption of k-means is wrong. Full covariances let it fit clusters that are stretched and tilted, of different sizes and different orientations, and the soft responsibilities let it be honest about points on a boundary instead of guessing. And because it's an actual probability model, it does things k-means can't: it gives you a density you can score new points against, it gives you a per-point confidence, and it gives you a principled way to choose the number of clusters. Add up the log-likelihood, penalize it for the number of parameters you spent, and you get the Bayesian information criterion — a score with a real minimum, so "how many clusters" stops being a judgment call about an elbow and becomes a number you can optimize.

The cost of that expressiveness is everything you'd expect from a bigger model. A full covariance in $D$ dimensions has $D(D+1)/2$ free parameters per component, so in high dimensions the mixture has a lot to estimate and can overfit, especially the component that grabs only a handful of points — its covariance can collapse toward a spike, sending the likelihood to infinity on a degenerate solution. The usual guard is a tiny ridge added to each covariance diagonal, which we do. EM also inherits the local-optimum problem from k-means and then some: the likelihood surface is bumpy, the answer depends on the initialization, and a bad start lands in a worse basin. The fix is the same — seed it well and restart it a few times.

The data

Three Gaussian blobs from scikit-learn's make_blobs, unequal sizes (110, 100, 90 points), then run through one shared linear map — a shear that stretches and tilts every blob into a long diagonal ellipse. That transform is deliberate: it's the canonical case where k-means fails and a full-covariance mixture succeeds, because the clusters are elongated in a direction that puts each one's far end closer to a neighbor's center than to its own. We commit the 300 points as data/blobs.csv along with the true blob id of each, which we treat as a sealed envelope: the mixture sees only the coordinates, and the ids come out at the end only to grade the clustering. The generator is code/make_data.py, seeded so the snapshot reproduces exactly.

Build it, one function at a time

Eight functions. The first is the shape term everything else calls; then the two EM steps; then the objective; then the seeding and the loop; then the two ways to read an answer out of a fitted model. Everything is pure NumPy.

Start with the multivariate Gaussian density. Given the points, a mean, and a covariance, we want the density at every point at once — no Python loop over the 300 rows:

def gaussian(X, mean, cov):
    """Multivariate Gaussian density N(x | mean, cov) for every row of X.

    X is (N, D), mean is (D,), cov is (D, D); the result is (N,), the density
    at each point. This is the shape term of the whole model — a full DxD
    covariance is what lets the level sets be tilted ellipses instead of the
    circles a single variance would give you.
    """
    D = X.shape[1]
    diff = X - mean                                   # (N, D)
    inv = np.linalg.inv(cov)                           # (D, D)
    # Mahalanobis distance squared for every row, without a Python loop.
    maha = np.einsum("ni,ij,nj->n", diff, inv, diff)   # (N,)
    norm = 1.0 / (np.power(2 * np.pi, D / 2) * np.sqrt(np.linalg.det(cov)))
    return norm * np.exp(-0.5 * maha)

The einsum computes the Mahalanobis distance for every row in one call, and the inv and det of the full covariance are what make the level sets tilted ellipses instead of circles. This one function is the entire difference from k-means; everything above it is bookkeeping.

The E-step stacks that density across components, weights each by its mixing weight, and normalizes each row so it sums to one:

def e_step(X, weights, means, covs):
    """E-step: the responsibility of every component for every point.

    Returns gamma, an (N, K) matrix where gamma[i, k] is the posterior
    probability that point i was drawn from component k, given the current
    parameters. Each row is a soft assignment that sums to 1 — this is the
    'soft' in soft clustering, and the whole reason a GMM differs from k-means.
    """
    K = len(weights)
    weighted = np.stack([weights[k] * gaussian(X, means[k], covs[k])
                         for k in range(K)], axis=1)   # (N, K): pi_k * N(x|k)
    total = weighted.sum(axis=1, keepdims=True)         # (N, 1): the mixture density
    return weighted / np.maximum(total, 1e-300)

The result gamma is the (N, K) matrix of responsibilities — the soft assignment. The 1e-300 floor keeps the division safe for a point so far from every component that every density underflows to zero. Given those responsibilities, the M-step re-fits every parameter as a weighted estimate:

def m_step(X, gamma, reg=1e-6):
    """M-step: re-estimate every parameter as a responsibility-weighted fit.

    Given the responsibilities, each component's new mean, covariance, and
    mixing weight is just a weighted version of the ordinary sample estimate,
    where every point contributes in proportion to how much it belongs to that
    component. `reg` adds a tiny ridge to each covariance diagonal so a
    component that grabs only a few points can't collapse to a singular matrix.
    """
    N, D = X.shape
    K = gamma.shape[1]
    Nk = gamma.sum(axis=0)                              # (K,): soft count per component
    weights = Nk / N                                    # mixing weights
    means = (gamma.T @ X) / Nk[:, None]                 # (K, D): weighted means
    covs = np.empty((K, D, D))
    for k in range(K):
        diff = X - means[k]                             # (N, D)
        covs[k] = (gamma[:, k, None] * diff).T @ diff / Nk[k]
        covs[k] += reg * np.eye(D)                       # keep it non-singular
    return weights, means, covs

The weighted mean is one matrix multiply; the weighted covariance is a short loop over the $K$ components, each one a responsibility-weighted scatter matrix straight from the math. The reg term is the ridge that keeps a small component from collapsing to a singular covariance. Now the number the loop climbs:

def log_likelihood(X, weights, means, covs):
    """Total log-likelihood of the data under the current mixture.

    Sum over points of the log of the mixture density. EM can only push this
    up or leave it flat, never down — it is the number the whole loop climbs,
    and the flat spot is where we stop.
    """
    K = len(weights)
    weighted = np.stack([weights[k] * gaussian(X, means[k], covs[k])
                         for k in range(K)], axis=1)   # (N, K)
    return float(np.log(np.maximum(weighted.sum(axis=1), 1e-300)).sum())

Same mixture density as the E-step, but instead of normalizing we take the log of the row sums and add them up. This can only rise across iterations; if it ever falls, there's a bug. Next the seeding, which matters here for the same reason it mattered in k-means:

def init_params(X, K, rng):
    """Seed the mixture: k-means++ means, one shared covariance, equal weights.

    Spreading the initial means with k-means++ (the same trick from the k-means
    chapter) keeps EM out of the bad local optima a careless random start falls
    into. Every component starts as one big round copy of the data's overall
    covariance, so the animation opens with three identical blobs that then peel
    apart and mold to their own cluster's shape.
    """
    # k-means++ seeding of the means
    first = int(rng.integers(len(X)))
    means = [X[first]]
    for _ in range(1, K):
        d2 = np.min([((X - m) ** 2).sum(axis=1) for m in means], axis=0)
        means.append(X[int(rng.choice(len(X), p=d2 / d2.sum()))])
    means = np.array(means, dtype=float)

    shared = np.cov(X.T)                                 # one covariance for all
    covs = np.stack([shared.copy() for _ in range(K)])
    weights = np.full(K, 1.0 / K)                        # equal mixing weights
    return weights, means, covs

We borrow k-means++ to spread the initial means, give every component one copy of the data's overall covariance, and start the weights equal. That's why the animation opens with three identical round blobs — they haven't found their own shapes yet. The loop ties it together: alternate E and M, recording a snapshot each pass, until the log-likelihood flattens:

def fit(X, K, rng, max_iter=100, tol=1e-4):
    """Run EM to convergence, recording a snapshot every iteration.

    Alternate the E-step and the M-step until the log-likelihood stops climbing
    by more than `tol`. Returns the final weights, means, covariances, and the
    full per-iteration history — each snapshot holds the parameters, the
    responsibilities, and the log-likelihood, which is exactly what the chapter
    replays frame by frame so you watch the ellipses form.
    """
    weights, means, covs = init_params(X, K, rng)
    gamma = e_step(X, weights, means, covs)
    ll = log_likelihood(X, weights, means, covs)
    history = [_snapshot(weights, means, covs, gamma, ll)]
    for _ in range(max_iter):
        weights, means, covs = m_step(X, gamma)         # re-estimate parameters
        gamma = e_step(X, weights, means, covs)         # recompute responsibilities
        new_ll = log_likelihood(X, weights, means, covs)
        history.append(_snapshot(weights, means, covs, gamma, new_ll))
        if new_ll - ll < tol:                           # likelihood stopped climbing
            ll = new_ll
            break
        ll = new_ll
    return weights, means, covs, ll, history

Each snapshot holds the full parameters, the responsibilities, and the log-likelihood — exactly what the animation replays frame by frame. Finally the two ways to read the fitted model. The hard label is the most responsible component:

def predict(X, weights, means, covs):
    """Hard label: the single most responsible component for each point.

    argmax over the responsibilities collapses the soft assignment back to one
    label per point, so a GMM can stand in for k-means when you need a firm
    clustering — but you threw away the confidence to get it.
    """
    return np.argmax(e_step(X, weights, means, covs), axis=1)

And the soft label is the whole responsibility vector — the thing k-means can never give you:

def predict_proba(X, weights, means, covs):
    """Soft label: the full responsibility vector for each point.

    This is what a GMM has that k-means never will — a distribution over
    components per point, so you can see the points on a boundary that belong
    honestly to two clusters at once.
    """
    return e_step(X, weights, means, covs)

Watch it work

This is the whole reason to stay in two dimensions. Below is a real run of the fit function above, one frame per EM iteration. Each component is drawn as two tilted ellipses — the 1σ and 2σ contours, computed from a real eigen-decomposition of that component's covariance, so the axes point along the true eigenvectors and their lengths are the square roots of the eigenvalues. Nothing here is a faked circle. The crosses are the means. The points are shaded by their soft responsibility: a point owned by one component takes its color, a point split between two comes out a blend, so the overlaps read as genuine mixtures. The lower panel tracks the log-likelihood climbing.

Press play. It opens at iteration 0 with three identical round contours sitting on the k-means++ seeds — one shared covariance, nobody committed yet. Then each frame does one E-step and one M-step, and you watch the ellipses peel apart, tilt into the diagonal of the data, and stretch to cover their own cluster while the points along the seams shift their blend. The log-likelihood shoots up in the first few frames and then crawls.

Watch the first three or four frames — that's where the ellipses find their angle, swinging from round to diagonal in a couple of steps. After that it's refinement: the contours nudge, the seam points settle their blends, and the log-likelihood climbs from about −1,081 at the seeded start to −515 at convergence, twenty-three iterations in, where the curve goes flat and the run stops. Reset and play it again — the seeding and everything else is seeded, so it molds the same three ellipses every time.

That repeatability is a property of this seed, not a promise. EM's likelihood surface is bumpy, and a different initialization can strand the ellipses in a worse arrangement — I hit exactly that on a couple of other seeds while building this, where the run locked two components onto one real cluster and split a third. That's the local-optimum caveat made concrete, and the reason the library restarts from several seeds by default.

The full implementation

The whole file, no library, top to bottom. This is exactly what the animation ran:

"""Gaussian mixture model with full covariances, fit by EM, from scratch.

k-means gave every point one hard label and drew every cluster as a round ball.
A Gaussian mixture keeps the idea of k centers but lets each cluster be a tilted,
stretched ellipse with soft edges: a point can belong 70% to one component and
30% to another. We fit it with expectation-maximization — the general method
lives in its own chapter; here it is specialized to the 2-D, full-covariance
mixture so we can watch the ellipses mold themselves to the data.

Pure NumPy. Every function below appears in the chapter one step at a time (the
`# region:` markers are what the book's include directives pull in).
"""

import numpy as np
import pandas as pd


# region: gaussian
def gaussian(X, mean, cov):
    """Multivariate Gaussian density N(x | mean, cov) for every row of X.

    X is (N, D), mean is (D,), cov is (D, D); the result is (N,), the density
    at each point. This is the shape term of the whole model — a full DxD
    covariance is what lets the level sets be tilted ellipses instead of the
    circles a single variance would give you.
    """
    D = X.shape[1]
    diff = X - mean                                   # (N, D)
    inv = np.linalg.inv(cov)                           # (D, D)
    # Mahalanobis distance squared for every row, without a Python loop.
    maha = np.einsum("ni,ij,nj->n", diff, inv, diff)   # (N,)
    norm = 1.0 / (np.power(2 * np.pi, D / 2) * np.sqrt(np.linalg.det(cov)))
    return norm * np.exp(-0.5 * maha)
# endregion


# region: e_step
def e_step(X, weights, means, covs):
    """E-step: the responsibility of every component for every point.

    Returns gamma, an (N, K) matrix where gamma[i, k] is the posterior
    probability that point i was drawn from component k, given the current
    parameters. Each row is a soft assignment that sums to 1 — this is the
    'soft' in soft clustering, and the whole reason a GMM differs from k-means.
    """
    K = len(weights)
    weighted = np.stack([weights[k] * gaussian(X, means[k], covs[k])
                         for k in range(K)], axis=1)   # (N, K): pi_k * N(x|k)
    total = weighted.sum(axis=1, keepdims=True)         # (N, 1): the mixture density
    return weighted / np.maximum(total, 1e-300)
# endregion


# region: m_step
def m_step(X, gamma, reg=1e-6):
    """M-step: re-estimate every parameter as a responsibility-weighted fit.

    Given the responsibilities, each component's new mean, covariance, and
    mixing weight is just a weighted version of the ordinary sample estimate,
    where every point contributes in proportion to how much it belongs to that
    component. `reg` adds a tiny ridge to each covariance diagonal so a
    component that grabs only a few points can't collapse to a singular matrix.
    """
    N, D = X.shape
    K = gamma.shape[1]
    Nk = gamma.sum(axis=0)                              # (K,): soft count per component
    weights = Nk / N                                    # mixing weights
    means = (gamma.T @ X) / Nk[:, None]                 # (K, D): weighted means
    covs = np.empty((K, D, D))
    for k in range(K):
        diff = X - means[k]                             # (N, D)
        covs[k] = (gamma[:, k, None] * diff).T @ diff / Nk[k]
        covs[k] += reg * np.eye(D)                       # keep it non-singular
    return weights, means, covs
# endregion


# region: log_likelihood
def log_likelihood(X, weights, means, covs):
    """Total log-likelihood of the data under the current mixture.

    Sum over points of the log of the mixture density. EM can only push this
    up or leave it flat, never down — it is the number the whole loop climbs,
    and the flat spot is where we stop.
    """
    K = len(weights)
    weighted = np.stack([weights[k] * gaussian(X, means[k], covs[k])
                         for k in range(K)], axis=1)   # (N, K)
    return float(np.log(np.maximum(weighted.sum(axis=1), 1e-300)).sum())
# endregion


# region: init_params
def init_params(X, K, rng):
    """Seed the mixture: k-means++ means, one shared covariance, equal weights.

    Spreading the initial means with k-means++ (the same trick from the k-means
    chapter) keeps EM out of the bad local optima a careless random start falls
    into. Every component starts as one big round copy of the data's overall
    covariance, so the animation opens with three identical blobs that then peel
    apart and mold to their own cluster's shape.
    """
    # k-means++ seeding of the means
    first = int(rng.integers(len(X)))
    means = [X[first]]
    for _ in range(1, K):
        d2 = np.min([((X - m) ** 2).sum(axis=1) for m in means], axis=0)
        means.append(X[int(rng.choice(len(X), p=d2 / d2.sum()))])
    means = np.array(means, dtype=float)

    shared = np.cov(X.T)                                 # one covariance for all
    covs = np.stack([shared.copy() for _ in range(K)])
    weights = np.full(K, 1.0 / K)                        # equal mixing weights
    return weights, means, covs
# endregion


# region: fit
def fit(X, K, rng, max_iter=100, tol=1e-4):
    """Run EM to convergence, recording a snapshot every iteration.

    Alternate the E-step and the M-step until the log-likelihood stops climbing
    by more than `tol`. Returns the final weights, means, covariances, and the
    full per-iteration history — each snapshot holds the parameters, the
    responsibilities, and the log-likelihood, which is exactly what the chapter
    replays frame by frame so you watch the ellipses form.
    """
    weights, means, covs = init_params(X, K, rng)
    gamma = e_step(X, weights, means, covs)
    ll = log_likelihood(X, weights, means, covs)
    history = [_snapshot(weights, means, covs, gamma, ll)]
    for _ in range(max_iter):
        weights, means, covs = m_step(X, gamma)         # re-estimate parameters
        gamma = e_step(X, weights, means, covs)         # recompute responsibilities
        new_ll = log_likelihood(X, weights, means, covs)
        history.append(_snapshot(weights, means, covs, gamma, new_ll))
        if new_ll - ll < tol:                           # likelihood stopped climbing
            ll = new_ll
            break
        ll = new_ll
    return weights, means, covs, ll, history
# endregion


# region: predict
def predict(X, weights, means, covs):
    """Hard label: the single most responsible component for each point.

    argmax over the responsibilities collapses the soft assignment back to one
    label per point, so a GMM can stand in for k-means when you need a firm
    clustering — but you threw away the confidence to get it.
    """
    return np.argmax(e_step(X, weights, means, covs), axis=1)
# endregion


# region: predict_proba
def predict_proba(X, weights, means, covs):
    """Soft label: the full responsibility vector for each point.

    This is what a GMM has that k-means never will — a distribution over
    components per point, so you can see the points on a boundary that belong
    honestly to two clusters at once.
    """
    return e_step(X, weights, means, covs)
# endregion


def _snapshot(weights, means, covs, gamma, ll):
    """One frame of the run: parameters, responsibilities, and log-likelihood."""
    return {
        "weights": weights.copy(),
        "means": means.copy(),
        "covs": covs.copy(),
        "gamma": gamma.copy(),
        "log_likelihood": ll,
    }


def load_data(path="../data/blobs.csv"):
    """Anisotropic blobs snapshot: 300 2-D points (x, y) plus the true blob id.

    The blob column is ground truth we NEVER fit on — the GMM sees only the
    coordinates. It exists so we can score the clustering afterwards.
    """
    return pd.read_csv(path)

The library version

Nobody hand-rolls a Gaussian mixture in production, and once you've built one you understand exactly what the library is doing. scikit-learn's GaussianMixture with covariance_type="full" is the same EM loop with full covariances, k-means initialization, several restarts so one unlucky seed can't decide the answer, and the same ridge guard against singular components. It optimizes the identical log-likelihood, so on this data it should land where we did:

def sklearn_gmm(X, k, seed=0):
    """Fit a full-covariance GMM with scikit-learn.

    Returns the hard labels, the component means, covariances, mixing weights,
    and the average per-sample log-likelihood (GaussianMixture.score) — the
    same things our fit() produces, so the two are directly comparable.
    """
    gm = GaussianMixture(n_components=k, covariance_type="full",
                         n_init=5, random_state=seed)
    labels = gm.fit_predict(X)
    return (labels, gm.means_, gm.covariances_, gm.weights_,
            float(gm.score(X)), gm)

The library also makes model selection cheap. Because the mixture is a real probability model, we can score each choice of $k$ with the Bayesian information criterion — the log-likelihood penalized by the parameter count — and pick the $k$ that minimizes it. Unlike the log-likelihood, which always improves as you add components, BIC has a genuine bottom:

def bic_sweep(X, ks, seed=0):
    """Fit a GMM for a range of k and record the Bayesian information criterion.

    Unlike log-likelihood, which always improves as you add components, BIC
    penalizes the parameter count — so it has a genuine minimum, and the k that
    minimizes it is the model-selection answer for 'how many clusters'.
    """
    out = []
    for k in ks:
        gm = GaussianMixture(n_components=k, covariance_type="full",
                             n_init=5, random_state=seed).fit(X)
        out.append((int(k), float(gm.bic(X))))
    return out

Fit the mixture for $k$ from 1 to 8 and plot the BIC of each:

BIC falls steeply from 1654.1 at k=1 to 1404.8 at k=2 to a minimum of 1126.7 at k=3, then turns and climbs — 1152.7 at k=4 and up from there, each extra component costing more in parameters than it earns in fit. The orange dot at k=3 is the answer, and it's the right one because we built the data with three blobs. This is the elbow question from k-means with a real objective under it: no eyeballing a bend, just the minimum of a curve.

Scratch versus library

Same data, same three components, both fit to convergence — our from-scratch EM against sklearn's GaussianMixture. The bars are the average per-point log-likelihood, the number both are maximizing:

They match to four decimals — both −1.7162. That's not rounding luck; both loops optimize the same likelihood on the same data and climb to the same optimum, and from a k-means++ start the basin is the good one. The two labelings agree perfectly too: adjusted Rand index 1.0 between our clusters and sklearn's, up to the arbitrary swap of which component is called 0, 1, or 2. And against the sealed envelope of true blob ids, both score a clean 1.0 — every point landed back in the group it was really drawn from.

Now the comparison this whole chapter is built around. Same 300 points, clustered two ways. On the left, k-means: hard labels, and a circle at each center reaching out to the cluster's spread. On the right, our GMM: hard labels from the same fit, each cluster wrapped in its 2σ tilted ellipse.

Look at the seams. k-means draws its circles and, because it judges every point by plain distance to a center, it cuts a clean straight line between the clusters — which runs diagonally across the real groups, mislabeling the elongated ends that poke past a neighbor's center. It scores an adjusted Rand index of 0.60 against the truth: a third of the points in the wrong group, held with total confidence. The GMM ellipses lie along the diagonal of each cluster, the Mahalanobis boundary bends around the shape, and the elongated ends stay with their own component. It scores 1.0. Same data, same k, same convergence guarantee — the only difference is that one model was allowed to draw an ellipse and the other wasn't.

Takeaways

A Gaussian mixture is what you reach for when k-means is the right instinct but the wrong shape. You still believe the data is a handful of blobs around centers; you just don't believe the blobs are round, equal, or cleanly separated. Full covariances buy you clusters that are stretched and tilted however the data demands, soft responsibilities buy you honesty at the boundaries, and because it's a real probability model you get things clustering alone never gives you: a density to score new points against, a per-point confidence, and BIC to pick the number of components without squinting at an elbow. On the data here that was the whole game — k-means at 0.60, the mixture at a perfect 1.0, decided entirely by whether a cluster was allowed to be an ellipse.

The bill comes due in two places, and both trace back to the same source: a full covariance is a lot of parameters. In high dimensions that's a lot to estimate and an easy way to overfit, and a component that grabs too few points can collapse its covariance to a spike and blow the likelihood up — which is why the ridge on the diagonal isn't optional. And EM lands in local optima, more of them than k-means, so the answer depends on the seed and you restart to be safe. When your clusters really are round and well-separated, k-means is faster and has less to go wrong; save the mixture for when the shape is the whole problem. It's the same trade the whole course keeps making — a more flexible model fits more, and asks you to be more careful about whether it fit or just memorized. The soft, shaped clustering here is also a doorway: the E-step-then-M-step loop, the responsibilities, the log-likelihood that only climbs, are the same moves behind hidden Markov models, topic models, and half of unsupervised learning. Build the mixture by hand once and you've built the pattern.