🎯 What You'll Learn
- Why K-Means' hard, equal-sized, spherical clusters break down on real data — and how a probabilistic model fixes it
- The Gaussian Mixture Model as a generative story: data is produced by K weighted Gaussians, and clustering means inferring which one produced each point
- The Expectation-Maximization (EM) algorithm: alternate between computing soft cluster "responsibilities" (E-step) and re-estimating each Gaussian's parameters (M-step)
- The four covariance types —
full,tied,diag,spherical— and the bias/flexibility trade-off each one makes - Choose the number of components with BIC/AIC instead of the elbow method, and know when GMM beats K-Means outright
1 From Hard to Soft Clustering: Why K-Means Isn't Enough
Lesson 28's limitations table ended with a specific promise: when clusters have different sizes, different shapes, or overlap, K-Means is the wrong tool. Look at what K-Means actually assumes, and the failure mode becomes obvious:
- Hard assignment. Every point belongs to exactly one cluster, full stop. A point sitting exactly between two natural groups still gets shoved into whichever centroid is a fraction closer — there is no way to say "this point is 60% cluster A, 40% cluster B."
- Implicitly spherical clusters. Because K-Means assigns by raw Euclidean distance to a centroid, its decision boundaries are always the perpendicular bisectors between centroids — clusters are forced to look like circles (in 2D) or spheres (in higher dimensions), never ellipses.
- Implicitly equal-sized, equal-spread clusters. A tight, dense cluster of 50 points and a loose, spread-out cluster of 500 points get treated identically — K-Means has no notion of "this cluster is naturally wider than that one."
A Gaussian Mixture Model (GMM) removes all three assumptions at once by replacing "nearest centroid" with "which probability distribution most likely generated this point?" Each cluster becomes its own Gaussian (bell-curve) distribution with its own mean, its own shape (covariance), and its own weight — and every point gets a soft probability of belonging to each cluster instead of a hard label.
import numpy as np
import matplotlib.pyplot as plt
from sklearn.datasets import make_blobs
from sklearn.cluster import KMeans
# Two clusters with very different spreads and sizes -- K-Means' nightmare
np.random.seed(42)
cluster_a = np.random.randn(300, 2) * [0.4, 0.4] + [0, 0] # tight, dense
cluster_b = np.random.randn(300, 2) * [2.2, 0.6] + [4, 0] # wide, elongated
X = np.vstack([cluster_a, cluster_b])
y_true = np.array([0] * 300 + [1] * 300)
km = KMeans(n_clusters=2, n_init=10, random_state=42)
km_labels = km.fit_predict(X)
# K-Means splits by perpendicular-bisector distance, not by cluster shape --
# it will steal points from the wide cluster's tail and hand them to the
# tight cluster, because "nearest centroid" ignores how spread out each
# cluster actually is.
from sklearn.metrics import adjusted_rand_score
print(f"K-Means ARI vs true clusters: {adjusted_rand_score(y_true, km_labels):.3f}")
# Meaningfully below 1.0 -- K-Means mis-assigns points near the tight
# cluster that actually belong to the wide one's near tail.K-Means asks "which centroid is closest?" GMM asks "which distribution was this point more likely sampled from?" That single reframing — geometry to probability — is what unlocks soft assignments, elliptical clusters, and clusters of different sizes, all for free.
2 The Gaussian Mixture Model: A Generative View of Clustering
A GMM assumes your data was generated by the following (imaginary, but useful) process, repeated for every data point:
- Pick one of K Gaussian components at random, with component k chosen with probability πk (the mixing weight — how common that cluster is overall). The weights sum to 1.
- Draw a point from that component's Gaussian distribution, N(μk, Σk) — its own mean vector μk and covariance matrix Σk.
The overall probability of observing a point x is the weighted sum over all K components:
p(x) = Σₖ πₖ · N(x | μₖ, Σₖ)
Fitting a GMM means finding the parameters — the K weights πk, means μk, and covariances Σk — that make the observed data most probable under this story. "Clustering" then becomes: for a given point, which component was most likely responsible for it?
import numpy as np
from scipy.stats import multivariate_normal
# A 2-component GMM's probability density is just a weighted sum of two
# multivariate normal densities -- nothing more exotic than that.
def gmm_pdf(x, weights, means, covariances):
density = 0.0
for pi_k, mu_k, sigma_k in zip(weights, means, covariances):
density += pi_k * multivariate_normal.pdf(x, mean=mu_k, cov=sigma_k)
return density
weights = [0.5, 0.5]
means = [np.array([0, 0]), np.array([4, 0])]
covariances = [np.eye(2) * 0.3, np.array([[2.0, 0], [0, 0.5]])]
point = np.array([2.0, 0.0]) # roughly halfway between the two components
print(f"p(x=[2,0]) under this GMM: {gmm_pdf(point, weights, means, covariances):.5f}")
# The "responsibility" of component k for a point is how much of that
# total density came from component k -- this is the soft cluster
# assignment, and it's the central quantity the EM algorithm computes.
for k, (pi_k, mu_k, sigma_k) in enumerate(zip(weights, means, covariances)):
contribution = pi_k * multivariate_normal.pdf(point, mean=mu_k, cov=sigma_k)
responsibility = contribution / gmm_pdf(point, weights, means, covariances)
print(f" Component {k}: responsibility = {responsibility:.3f}")If you force every component's covariance to be spherical and identical (Σk = σ²I for all k) and take the limit as σ→0, soft responsibilities collapse into hard 0/1 assignments and GMM's parameter updates become exactly K-Means' centroid updates. K-Means isn't a different algorithm from GMM — it's GMM with the flexibility deliberately switched off.
3 The EM Algorithm: Expectation and Maximization
There's a chicken-and-egg problem in fitting a GMM: if you knew which component generated each point, estimating each Gaussian's mean and covariance would be trivial (just compute the mean/covariance of its assigned points). And if you knew each Gaussian's parameters, computing each point's responsibilities would be trivial (Section 2's formula). But you know neither at the start.
Expectation-Maximization (EM) breaks the deadlock by alternating between the two, exactly like K-Means alternates between assignment and update — just with soft probabilities instead of hard labels:
- Initialize the K components' weights, means, and covariances (commonly by running K-Means first and using its clusters as a starting point).
- E-step (Expectation): With the current parameters fixed, compute every point's responsibility for every component — "given what I currently believe about the clusters, how likely is each cluster to have generated this point?"
- M-step (Maximization): With the responsibilities fixed, re-estimate each component's weight, mean, and covariance as a responsibility-weighted average over all points, instead of a hard-membership average over assigned points only.
- Repeat E and M until the total log-likelihood of the data stops improving (convergence).
Just like K-Means, each EM iteration is guaranteed to never decrease the data's log-likelihood, so the algorithm always converges. But the likelihood surface has many local maxima, and a bad random initialization can strand you at a mediocre one. This is why sklearn.mixture.GaussianMixture exposes n_init — run the whole EM procedure several times from different starts and keep the run with the highest final log-likelihood, the same defense K-Means uses.
Watch It Run: EM Step by Step
The simulator below fits a 2-component GMM to the same "one tight cluster, one wide cluster" dataset from Section 1 — exactly the case where K-Means struggled. Each point's color is a blend of the two cluster colors, weighted by its current responsibilities (a point that's 70% likely to be from cluster A and 30% likely from cluster B is rendered as a 70/30 color mix). The ellipses show each component's current covariance — one standard deviation out. Click Step to alternate E-steps and M-steps, and watch the ellipses stretch to match the true cluster shapes while K-Means-style circles never could.
Iteration 0 — components initialized from a K-Means warm start. Click Step to run the E-step and M-step.
4 E-Step: Computing Responsibilities
The responsibility γ(znk) — "how responsible is component k for point n" — is computed with exactly the Bayes'-theorem-flavored formula from Section 2, applied to every point against every component simultaneously:
γ(zₙₖ) = πₖ N(xₙ | μₖ, Σₖ) / Σⱼ πⱼ N(xₙ | μⱼ, Σⱼ)
The numerator is "how likely is component k, weighted by how common it is, to have produced this exact point." The denominator normalizes across all K components so that a single point's responsibilities always sum to 1 — every point's probability mass is fully accounted for, split across the components in proportion to how well each one explains it.
import numpy as np
from scipy.stats import multivariate_normal
def e_step(X, weights, means, covariances):
"""Compute the (n_samples, n_components) responsibility matrix."""
n_samples, K = X.shape[0], len(weights)
resp = np.zeros((n_samples, K))
for k in range(K):
resp[:, k] = weights[k] * multivariate_normal.pdf(
X, mean=means[k], cov=covariances[k])
# Normalize each row to sum to 1 (Bayes' theorem denominator)
resp = resp / resp.sum(axis=1, keepdims=True)
return resp
np.random.seed(42)
X_demo = np.vstack([
np.random.randn(5, 2) * 0.3 + [0, 0],
np.random.randn(5, 2) * 0.3 + [4, 0],
])
weights = [0.5, 0.5]
means = [np.array([0.2, 0.1]), np.array([3.8, -0.1])]
covariances = [np.eye(2) * 0.4, np.eye(2) * 0.4]
resp = e_step(X_demo, weights, means, covariances)
print("Responsibilities (rows sum to 1):")
print(resp.round(3))
# Points near [0,0] get responsibility ~[0.99, 0.01] for component 0;
# points near [4,0] get roughly the mirror image.The responsibility matrix is the clustering result — a full probability distribution over components for every point, not a single label. GaussianMixture.predict_proba(X) returns exactly this matrix; .predict(X) just takes the argmax of it for you when you need a hard label.
5 M-Step: Updating the Parameters
The M-step re-estimates each component's parameters as a weighted version of ordinary maximum-likelihood estimates — every point contributes to every component's update, weighted by that point's responsibility for that component. Define Nk = Σn γ(znk) — the "effective number of points" softly assigned to component k:
import numpy as np
def m_step(X, resp):
"""Re-estimate weights, means, and covariances from responsibilities."""
n_samples, K = resp.shape
N_k = resp.sum(axis=0) # effective count per component
weights = N_k / n_samples # new mixing weights: pi_k = N_k / n
means = (resp.T @ X) / N_k[:, None] # responsibility-weighted mean
covariances = []
for k in range(K):
diff = X - means[k] # (n_samples, n_features)
# Responsibility-weighted outer product, summed over points
weighted_cov = (resp[:, k:k+1] * diff).T @ diff / N_k[k]
covariances.append(weighted_cov)
return weights, means, covariances
# Continuing the E-step example: feed its responsibilities back through M
from scipy.stats import multivariate_normal
def e_step(X, weights, means, covariances):
n_samples, K = X.shape[0], len(weights)
resp = np.zeros((n_samples, K))
for k in range(K):
resp[:, k] = weights[k] * multivariate_normal.pdf(X, mean=means[k], cov=covariances[k])
return resp / resp.sum(axis=1, keepdims=True)
np.random.seed(42)
X_demo = np.vstack([np.random.randn(50, 2) * 0.4 + [0, 0], np.random.randn(50, 2) * 0.4 + [4, 0]])
weights, means, covariances = [0.5, 0.5], [np.array([1.0, 0.0]), np.array([3.0, 0.0])], [np.eye(2), np.eye(2)]
for iteration in range(10):
resp = e_step(X_demo, weights, means, covariances)
weights, means, covariances = m_step(X_demo, resp)
print("Converged weights:", np.round(weights, 3))
print("Converged means:\n", np.round(means, 3))
# Should land close to pi=[0.5, 0.5], means near [0,0] and [4,0]Notice the parallel to K-Means' update step: K-Means recomputes each centroid as the mean of its hard-assigned points; GMM recomputes each component's mean (and now also its covariance) as the mean of all points, weighted by responsibility. Points far from a component barely move its estimate; points close to it dominate. The M-step is doing weighted maximum likelihood estimation, one component at a time.
6 Covariance Types: full, diag, tied, spherical
The covariance matrix Σk is what lets a component be an ellipse instead of a circle — but a full covariance matrix has a lot of parameters (for d features, d(d+1)/2 per component), which can overfit on small datasets or high dimensions. scikit-learn's GaussianMixture exposes a covariance_type argument that trades flexibility for parameter count:
| covariance_type | Shape it allows | Parameters per component |
|---|---|---|
full (default) |
Any ellipse — any size, shape, and orientation, independent per component | d(d+1)/2 — most flexible, most prone to overfitting |
tied |
Any ellipse, but every component shares the same shape and orientation | d(d+1)/2 total (shared across all K) |
diag |
Axis-aligned ellipse (no tilt) — each feature has its own spread, independently per component | d per component |
spherical |
A circle/sphere — one radius per component, no elongation at all | 1 per component (closest thing to K-Means) |
import numpy as np
from sklearn.mixture import GaussianMixture
from sklearn.datasets import make_blobs
X, y_true = make_blobs(n_samples=600, centers=3, cluster_std=[0.5, 1.5, 0.8], random_state=42)
for cov_type in ['full', 'tied', 'diag', 'spherical']:
gmm = GaussianMixture(n_components=3, covariance_type=cov_type, n_init=5, random_state=42)
gmm.fit(X)
print(f"{cov_type:10s} log-likelihood: {gmm.score(X) * len(X):10.1f} "
f"BIC: {gmm.bic(X):8.1f} n_params: {gmm._n_parameters()}")
# 'full' usually gets the best (lowest) BIC when clusters genuinely differ in
# shape, like this dataset -- but on small or high-dimensional data,
# 'diag' or 'tied' can win by avoiding overfitting the extra parameters.Start with covariance_type='full' — it's the most expressive and sklearn's default. If you have very high-dimensional data or suspect overfitting (fit looks great on training data, unstable across re-runs), compare BIC across all four covariance types (Section 7) and let the data pick the right trade-off rather than guessing.
7 Choosing K: BIC and AIC
K-Means' elbow method doesn't transfer to GMM, because a GMM's log-likelihood — unlike WCSS — doesn't monotonically improve in a way that produces a clean elbow; it can also be pushed arbitrarily high by overfitting (imagine one tiny-variance Gaussian centered exactly on each point). Instead, GMM model selection uses information criteria that explicitly penalize complexity:
- BIC (Bayesian Information Criterion): −2·log-likelihood + k·log(n), where k is the number of parameters and n is the number of samples. The log(n) penalty makes BIC increasingly strict about extra parameters as your dataset grows.
- AIC (Akaike Information Criterion): −2·log-likelihood + 2·k. A softer, constant penalty per parameter regardless of dataset size — tends to prefer slightly larger, more flexible models than BIC.
For both, lower is better — you're looking for the number of components that best balances how well the model fits the data against how many parameters it took to get there.
import numpy as np
import matplotlib.pyplot as plt
from sklearn.mixture import GaussianMixture
from sklearn.datasets import make_blobs
X, _ = make_blobs(n_samples=600, centers=4, cluster_std=0.9, random_state=42)
n_components_range = range(1, 11)
bic_scores, aic_scores = [], []
for n in n_components_range:
gmm = GaussianMixture(n_components=n, covariance_type='full', n_init=5, random_state=42)
gmm.fit(X)
bic_scores.append(gmm.bic(X))
aic_scores.append(gmm.aic(X))
print(f"K={n:2d} BIC={gmm.bic(X):9.1f} AIC={gmm.aic(X):9.1f}")
best_k_bic = list(n_components_range)[np.argmin(bic_scores)]
best_k_aic = list(n_components_range)[np.argmin(aic_scores)]
print(f"\nBest K by BIC: {best_k_bic} (true K=4)")
print(f"Best K by AIC: {best_k_aic}")
plt.figure(figsize=(9, 5))
plt.plot(n_components_range, bic_scores, 'o-', label='BIC', linewidth=2)
plt.plot(n_components_range, aic_scores, 's-', label='AIC', linewidth=2)
plt.axvline(x=best_k_bic, color='red', linestyle='--', alpha=0.6, label=f'BIC minimum (K={best_k_bic})')
plt.xlabel('Number of Components (K)')
plt.ylabel('Information Criterion (lower = better)')
plt.title('Choosing K for a GMM with BIC / AIC')
plt.legend()
plt.tight_layout()
plt.savefig('gmm_bic_aic.png', dpi=150)If your true data-generating process has 4 clusters but two of them overlap heavily, BIC may well prefer K=3 because a 3-component model explains the data almost as well with fewer parameters. Treat the BIC-optimal K as a strong statistical recommendation, then sanity-check it against domain knowledge before locking it in.
8 GMM vs K-Means: When to Use Which
import numpy as np
from sklearn.cluster import KMeans
from sklearn.mixture import GaussianMixture
from sklearn.datasets import make_blobs
from sklearn.metrics import adjusted_rand_score
import time
# Recreate the "tight cluster + wide cluster" scenario from Section 1
np.random.seed(42)
cluster_a = np.random.randn(300, 2) * [0.4, 0.4] + [0, 0]
cluster_b = np.random.randn(300, 2) * [2.2, 0.6] + [4, 0]
X = np.vstack([cluster_a, cluster_b])
y_true = np.array([0] * 300 + [1] * 300)
t0 = time.time()
km_labels = KMeans(n_clusters=2, n_init=10, random_state=42).fit_predict(X)
km_time = time.time() - t0
t0 = time.time()
gmm_labels = GaussianMixture(n_components=2, covariance_type='full', n_init=5, random_state=42).fit_predict(X)
gmm_time = time.time() - t0
print(f"K-Means ARI: {adjusted_rand_score(y_true, km_labels):.3f} ({km_time*1000:.1f} ms)")
print(f"GMM ARI: {adjusted_rand_score(y_true, gmm_labels):.3f} ({gmm_time*1000:.1f} ms)")
# GMM's ARI is noticeably higher here -- it correctly models the wide
# cluster's larger footprint instead of splitting it against a same-size
# assumption. The cost: GMM is doing more work per iteration (covariance
# estimation and matrix inversion), so it's meaningfully slower per fit.| Use K-Means when... | Use GMM when... |
|---|---|
| Clusters are roughly spherical and similar in size (common after good feature scaling) | Clusters are elongated, differently sized, or overlap |
You need speed on very large datasets (millions of points) — use MiniBatchKMeans |
You need calibrated probabilities, not just labels — e.g. "60% confident this transaction is fraud-cluster" |
| You just need a quick, interpretable baseline for segmentation | You want a generative model — e.g. to sample new synthetic points, or score how "typical" a new point is (anomaly detection, Lesson 33) |
| Data is high-dimensional and covariance estimation would be unstable | You want soft, principled model comparison via BIC/AIC instead of a heuristic elbow |
In practice, a very common workflow is to run K-Means first — it's fast and gives a decent starting guess — and use its cluster assignments to initialize a GMM's means, letting EM refine the boundaries into whatever shapes the data actually has. sklearn.mixture.GaussianMixture does exactly this by default (init_params='kmeans').
Real-World Spotlight: Soft Fraud Segmentation with a GMM
A payments company wants to segment transactions into behavioral clusters for downstream fraud modeling — but a hard K-Means label throws away useful information: a transaction sitting near the boundary between "normal spending" and "unusual burst" is exactly the kind of ambiguous case a fraud analyst wants flagged for review, not silently assigned to whichever cluster happens to win.
import numpy as np
import pandas as pd
from sklearn.mixture import GaussianMixture
from sklearn.preprocessing import StandardScaler
np.random.seed(42)
n = 3000
# Simulate transaction features: amount, velocity (txns/hour), merchant_risk_score
df = pd.DataFrame({
'amount': np.concatenate([
np.random.lognormal(3.5, 0.6, 2600), # typical spending
np.random.lognormal(6.0, 0.9, 400), # unusual bursts
]),
'velocity': np.concatenate([
np.random.poisson(1.2, 2600),
np.random.poisson(6.0, 400),
]),
'merchant_risk_score': np.concatenate([
np.random.beta(2, 8, 2600),
np.random.beta(6, 3, 400),
]),
})
features = ['amount', 'velocity', 'merchant_risk_score']
X_scaled = StandardScaler().fit_transform(df[features])
# ── Select K with BIC ──
bic_scores = []
for k in range(1, 6):
gmm = GaussianMixture(n_components=k, covariance_type='full', n_init=5, random_state=42)
gmm.fit(X_scaled)
bic_scores.append(gmm.bic(X_scaled))
best_k = np.argmin(bic_scores) + 1
print(f"Best K by BIC: {best_k}")
# ── Fit final GMM and inspect soft assignments ──
gmm = GaussianMixture(n_components=best_k, covariance_type='full', n_init=10, random_state=42)
df['cluster'] = gmm.fit_predict(X_scaled)
probs = gmm.predict_proba(X_scaled)
# The key advantage over K-Means: flag ambiguous transactions for review
# instead of silently hard-assigning them.
df['max_prob'] = probs.max(axis=1)
ambiguous = df[df['max_prob'] < 0.65]
print(f"\nTransactions with no dominant cluster (max prob < 0.65): {len(ambiguous)} "
f"({len(ambiguous) / len(df):.1%})")
print("These are prime candidates for manual fraud-analyst review --")
print("a hard-labeling algorithm like K-Means would have silently picked one.")
# ── Profile the clusters by their means (original scale) ──
scaler = StandardScaler().fit(df[features])
cluster_means = pd.DataFrame(
scaler.inverse_transform(gmm.means_), columns=features
).round(2)
cluster_means['weight'] = gmm.weights_.round(3)
print("\nCluster profiles:")
print(cluster_means)The GMM's soft probabilities do real work here: instead of a binary "fraud cluster / not fraud cluster" label, every transaction gets a full probability distribution over behavioral segments, and the ones with no dominant segment become a natural, principled review queue — a distinction K-Means' hard labels simply cannot express.
✍️ Practice Exercises
- Generate two overlapping 2D Gaussian clusters with
make_blobsusing differentcluster_stdvalues for each. Fit bothKMeansandGaussianMixturewith K=2. Compare theiradjusted_rand_scoreagainst the true labels — by how much does GMM win, and why? - On the Iris dataset, fit a GMM with
covariance_type='full'and K=3. Usepredict_probato find the 5 flowers with the lowest maximum responsibility (the most "ambiguous" cases). Are they near a species boundary? - Fit GMMs with K from 1 to 8 on a dataset of your choice, plotting BIC and AIC on the same axes. Do they agree on the best K? If not, which one matches your domain intuition better?
- Compare all four
covariance_typesettings on a small (n=50), high-dimensional (10+ feature) synthetic dataset. Which one achieves the best BIC, and why doesfulltend to struggle here specifically?
▶ Hints
from sklearn.datasets import load_iris
from sklearn.mixture import GaussianMixture
import numpy as np
iris = load_iris()
gmm = GaussianMixture(n_components=3, covariance_type='full', n_init=10, random_state=42)
gmm.fit(iris.data)
probs = gmm.predict_proba(iris.data)
max_probs = probs.max(axis=1)
most_ambiguous = np.argsort(max_probs)[:5]
print(iris.data[most_ambiguous])
print(iris.target[most_ambiguous]) # check which species these belong to📚 Primary Source for This Lesson
scikit-learn: Gaussian Mixture Models
The official guide covers all four covariance types, model selection with BIC/AIC, and the Bayesian (Dirichlet Process) variant. For the original EM algorithm, see Dempster, Laird & Rubin (1977) "Maximum Likelihood from Incomplete Data via the EM Algorithm" (Journal of the Royal Statistical Society) — the paper that gave the general E-step/M-step framework its name, of which GMM fitting is the most common special case.