Real-world heights aren't one bell curve — they're two (men and women) overlapping. A single Gaussian gets the mean right and the variance terrible. GMMs fix this. They're also the foundation of every soft-clustering algorithm you'll ever use.
Learning Objectives
After this lesson, you will be able to:
Model data as a mixture of K Gaussian components — soft clusters where each point belongs partially to several clusters, weighted by responsibilities, instead of getting a single hard label like in K-Means
Run the Expectation-Maximization (EM) algorithm by hand and in code — alternating between an E-step that computes responsibilities and an M-step that updates means, covariances, and mixing weights to maximize the data likelihood
Pick the right covariance type (spherical, diagonal, tied, full) and the right number of components K using BIC or AIC, instead of guessing
Use GMM beyond clustering — for density estimation, anomaly detection (low-likelihood points are outliers), and generative sampling — and know when to switch to DBSCAN or HDBSCAN because the clusters are non-elliptical
Don't worry if the math feels heavier than K-Means at first — once you see the E-step and M-step alternating, EM clicks fast and you have a tool that works on far more than just clustering.
Playground:Clustering → — the same playground supports GMM mode; flip between K-Means hard assignments and GMM soft probabilities to see the difference live.
Before going further, get your hands on the underlying primitive — a single Gaussian and a sum-of-Gaussians — using the distribution viz to see how mixing weight, mean, and variance shape the resulting density curve.
Single Gaussian vs mixture of Gaussians — adjust the components and watch the density changeInteractive
Loading visualization...
Try this: Start with a single Gaussian and shift its mean / stretch its variance. Then enable a second component with a different mean — watch the combined density form two peaks. That's exactly what a 2-component GMM is doing: tuning the means, variances, and mixing weights of two Gaussians until their weighted sum best explains the data.
A Gaussian Mixture Model assumes your data was generated by the following recipe: there are K hidden Gaussian distributions in the population, with mixing weights π_1 ... π_K that sum to 1. To produce a data point, nature flips a weighted K-sided coin to pick component k, then samples a point from that component's Gaussian N(μ_k, Σ_k). You see only the points -- not which component produced them. Your job is to recover the K means, K covariances, and K mixing weights from the data alone.
In K-Means, point i belongs to cluster k or it does not -- a 0 or a 1. In GMM, each point has a responsibility vector -- a probability distribution over which cluster generated it.
γik=∑j=1KπjN(xi∣μj,Σj)πkN(xi∣μk,Σk)
The responsibilities are the bridge between hard and soft clustering. If you set the largest γ to 1 and the rest to 0, you recover hard cluster assignments -- this is why K-Means is sometimes called the "hard-EM" limit of GMM with spherical equal-variance covariance.
Given the current parameter estimates , compute γ_ik for every point i and every cluster k using the formula above. This step is "expectation" because it computes the expected component membership of each point under the current parameters.
Given the current responsibilities, update the parameters to maximize the expected complete-data likelihood. Each cluster's parameters are weighted averages where the weights are the responsibilities.
After each EM round, the data log-likelihood is guaranteed not to decrease. Either you climb or you plateau -- you never slide backward. This is the key theoretical property of EM. The catch: you converge to a local maximum of the likelihood, not necessarily the global one. The fix: run EM from multiple random starts (sklearn defaults to n_init=1 -- you should usually set it to 5 or 10) and keep the best.
The "log-likelihood never decreases" claim isn't a coincidence — it falls out of a more general object called the Evidence Lower BOund (ELBO). Once you see EM through the ELBO lens, monotonic improvement is obvious, and you also have the bridge to variational inference, VAEs, and modern Bayesian deep learning.
The trouble is that ln p(X | θ) involves a sum inside a log (marginalizing out the latent cluster assignments Z), which is intractable to maximize directly. The trick is to introduce an arbitrary distribution q(Z) over the latents and apply Jensen's inequality:
This single inequality is the foundation of EM. Both steps now have crisp interpretations:
E-step (tighten the bound). With θ held fixed, the choice of q(Z) that makes the ELBO largest — and therefore closes the gap completely — is the posterior q*(Z) = p(Z | X, θ). That's exactly the responsibility computation: γ_ik = P(component k | x_i, current θ). After the E-step, ELBO(q*, θ) = ln p(X | θ).
M-step (push the bound up). Now hold q fixed and maximize the ELBO over θ. Since the entropy term doesn't depend on θ, this collapses to maximizing the expected complete-data log-likelihood — the weighted MLE we wrote down for μ_k, Σ_k, π_k. Because the ELBO is a lower bound on ln p(X | θ) and the E-step just made it equal to the log-likelihood, any M-step improvement directly improves the log-likelihood.
E-step closes the gaplnp(X∣θ(t))=ELBO(q(t+1),θ(t))≤M-step pushes upELBO(q(t+1),θ(t+1))≤lnp(X∣θ(t+1))
What Do You Think?
K-Means and GMM are both run on the same 2D dataset. K-Means converges to silhouette=0.32, GMM converges to a much higher data log-likelihood. What does the data most likely look like?
The answer is option B. K-Means assumes spherical, equal-size clusters and gets confused when clusters are elongated or unequal. GMM with full covariance fits ellipsoids of any shape and any size, so it dominates K-Means on this kind of data. For interlocking half-moons (option C), neither GMM nor K-Means works -- you need DBSCAN or spectral clustering.
sklearn.mixture.GaussianMixture(covariance_type=...) exposes four choices, and picking the wrong one is the second-most-common GMM mistake (after n_init=1).
Type
Σ_k shape
Free parameters
Best for
spherical
σ_k² · I
K
Tiny data; almost-K-Means
diagonal
diag(σ²_, ..., σ²_)
K·d
Independent features per cluster
tied
shared Σ across clusters
d·(d+1)/2
Linear discriminant flavor; clusters parallel
full
unique Σ_k per cluster
K·d·(d+1)/2
Generic case; default if you have data
Default to full when N is at least 50 × the number of free parameters. Drop to diagonal when you have many features and not much data. tied is essentially Linear Discriminant Analysis -- useful when you believe all clusters have the same shape.
Unlike K-Means with elbow plots, GMM gives you a real likelihood number, which lets you use principled model-selection criteria.
BIC=−2lnL^+plnN
In practice, plot BIC and AIC across a grid of K values. The knee of the BIC curve is your best bet -- one component above and below the knee should give noticeably worse fit. If BIC keeps decreasing with K, your data may not actually be a mixture of Gaussians (try DBSCAN). If BIC rises immediately, you do not have multiple components.
#Worked Example: One Round of EM by Hand (1D, K=2)
Three data points: x = [1, 2, 8]. Initial guesses: π = (0.5, 0.5), μ₁ = 0, μ₂ = 10, σ₁ = σ₂ = 2.
E-step. For each point, compute γ_ik = π_k · N(x | μ_k, σ_k²) / Σ_j …
σ₁² = (1·(1−1.5)² + 1·(2−1.5)² + 0·…)/2 = 0.25, so σ₁ = 0.5. σ₂² = 1·(8−8)²/1 = 0, but in real implementations a tiny ridge reg_covar=1e-6 is added so σ₂ doesn't collapse.
One iteration moved μ₁ from 0 → 1.5 and μ₂ from 10 → 8 — both much closer to where the data actually lives. The next E-step recomputes responsibilities under the new parameters, and so on until log-likelihood plateaus.
Run a real GaussianMixture on a 3-blob dataset and watch BIC pick the right K while soft probabilities for boundary points reveal the cluster ambiguity that K-Means would have hidden.
Loading visualization...
Quick check
In the GMM result above, a point sits exactly between two blob centers and gets soft probabilities (0.51, 0.49, 0.00). What does this tell you, and how would K-Means handle the same point?
Quick check
During the M-step of EM for a GMM, you update the mean of component k. The formula is the responsibility-weighted average of all data points. Why are points multiplied by their responsibility γ_ik instead of just summed?
Try it: see how EM converges from random initial parameters to fitted Gaussian componentsInteractive
Loading visualization...
EM is famously fast in low dimensions -- typically 10-50 iterations to convergence on well-separated data, and 100-200 on heavy overlap. The runtime per iteration is O(N · K · d²) for full covariance, so it scales linearly in dataset size and quadratically in feature dimension.
GMM gives you p(x) for any new point x. Plot p(x) over a grid and you have a full density estimate of your data -- useful for sampling, simulation, and anomaly detection.
Train a GMM on normal data only. For a new point x, the log-likelihood ln p(x) under the trained GMM tells you how "normal" x is. Threshold the log-likelihood and you have an unsupervised anomaly detector. This is the GMM-based approach to outlier detection -- far more principled than the Z-score for multivariate, non-Gaussian-marginal data.
Sampling from a GMM is one line: pick a component k with probability π_k, then sample x ~ N(μ_k, Σ_k). This is the simplest generative model in the classical-ML toolkit, and it underlies many of the synthetic-data techniques used before GANs and diffusion models took over for images.
from sklearn.mixture import GaussianMixture
from sklearn.datasets import make_blobs
import numpy as np
X, y_true = make_blobs(n_samples=500, centers=3, cluster_std=1.5, random_state=0)
# COMMON BUG: n_init=1 (the default) → bad local optima on most runs
# FIX: always use n_init=5 or more
gmm = GaussianMixture(n_components=3, covariance_type='full', n_init=10, random_state=0)
gmm.fit(X)
print(f"Log-likelihood: {gmm.score(X) * len(X):.1f}")
print(f"BIC: {gmm.bic(X):.1f}")
print(f"Mixing weights: {gmm.weights_.round(3)}")
# Expected: mixing weights near [0.33, 0.33, 0.33], BIC noticeably lower than K=2 or K=4
For 3 well-separated balanced blobs, mixing weights should land near (0.33, 0.33, 0.33) and BIC should be minimized at K=3 across a sweep K=1..6.
#GMM vs. K-Means vs. DBSCAN vs. HDBSCAN: When Each Wins
Situation
Reach for
Spherical, equal-size clusters, hard labels OK
K-Means (faster, simpler)
Elliptical/rotated clusters, want soft probabilities or p(x)
GMM (covariance_type='full')
Cluster shapes you can't predict, noise/outliers present
DBSCAN
DBSCAN but density varies across regions
HDBSCAN
Don't know K, want to discover it
BayesianGaussianMixture (Dirichlet process) or HDBSCAN
Need density estimation for anomaly detection
GMM (negative log-likelihood = anomaly score)
Need to sample synthetic data from clusters
GMM's superpower over K-Means: it gives you p(x) for free, which means density estimation, anomaly detection, and generative sampling all come bundled — three jobs for the price of one fit.
Quick check
A dataset has three elongated, rotated, partially overlapping clusters with different sizes. Compare GMM (covariance_type='full') vs K-Means on this data.
GMM is soft clustering with a generative model behind it. Every point gets a probability distribution over components, every component is a full Gaussian with mean and covariance, and the model gives you p(x) for any new point so you can do density estimation, anomaly detection, and sampling -- not just clustering
EM alternates a chicken-and-egg loop until convergence. E-step computes responsibilities given current parameters using Bayes' rule, M-step updates parameters as responsibility-weighted averages, and the data log-likelihood never decreases between iterations
The covariance type is a major hyperparameter. full is the default flexible choice, diagonal for high-d when you trust feature independence within clusters, tied for LDA-flavored shared shape, spherical only for tiny data; mismatching the covariance type to your data shape cripples the fit
Always set n_init at least 5 and pick K with BIC. Single-init EM lands on bad local optima frequently; BIC's penalty on parameters protects you from over-fitting K, AIC is more permissive; if BIC keeps falling forever, your data probably is not a Gaussian mixture
Switch to DBSCAN or HDBSCAN when clusters are non-convex. GMM forces ellipsoids; rings, half-moons, and manifolds need density-based or spectral methods, no amount of K tuning rescues GMM on those shapes
What is the main difference between K-Means and a Gaussian Mixture Model?
GMMs gave you soft clusters and a real density model. Next up: Anomaly Detection — using density estimates, isolation forests, and one-class SVMs to find the rare points that don't fit anywhere.