Gaussian mixture models

A marketing team has 100,000 customers and two useful columns: weekly app visits and average order value. A customer who visits twice a week is probably an occasional buyer. A customer who visits ten times is probably a regular.

The awkward customers sit between those stories. Someone with seven visits might be a regular having a quiet month, or an occasional buyer during a promotion. K-means must choose one label. That choice is convenient, but it can pretend to know more than the data does.

A Gaussian mixture model, or GMM, models the data as several overlapping bell-shaped distributions. It gives each observation a posterior membership probability: the model’s probability that the observation came from each component, rather than a forced hard label. It can also learn clusters that are stretched, tilted, or different in size.

The word mixture is literal:

1. Choose a hidden group according to its mixture weight.
2. Generate a data point from that group’s Gaussian distribution.

The hidden group is a component. It is not automatically a real-world customer type; it is one Gaussian distribution that helps explain the observed data.

Soft assignment, with arithmetic

Suppose we temporarily use only weekly visits. The model has two components:

• Component A has weight 0.65, mean 4 visits, and standard deviation 2.
• Component B has weight 0.35, mean 10 visits, and standard deviation 2.

A mixture weight is the proportion of the model assigned to a component. The weights add to 1; the mean is the centre, and the standard deviation describes typical spread.

Consider a customer with 8 visits. The Gaussian density is about 0.027 under A and 0.121 under B. B is a better explanation because 8 is one standard deviation below B’s mean but two above A’s.

Include the mixture weights:

• A contributes 0.65 × 0.027 = 0.0176.
• B contributes 0.35 × 0.121 = 0.0424.
• The total mixture density is 0.0600.

Normalising those contributions gives:

• A: 0.0176 / 0.0600 = 0.293, or about 29%.
• B: 0.0424 / 0.0600 = 0.707, or about 71%.

So predict_proba would report roughly 29% for A and 71% for B, subject to the fitted parameters. K-means would report one cluster ID and discard the uncertainty.

This probability answers a conditional question: given this fitted model and observation, which component is the more plausible source?

The overlap supports a cautious rule such as “send the regular-buyer offer only when B exceeds 0.8.” Customers around 0.5 can be held out instead of quietly mislabelled.

What the model is fitting

In several dimensions, a Gaussian becomes an ellipse-shaped density. Its centre is a mean vector, and its covariance matrix records feature spread and whether features move together.

For component k, the parameters are:

• μ_k: mean vector
• Σ_k: covariance matrix
• π_k: mixture weight

The mixture density is:

p(x) = Σ π_k N(x | μ_k, Σ_k)

A nearby component contributes more density; so does one with a larger weight. The membership probability, or responsibility, is:

r_ik = π_k N(x_i | μ_k, Σ_k) / Σ_j π_j N(x_i | μ_j, Σ_j)

The denominator makes responsibilities add to 1 for each observation. This is the mechanism behind soft clustering, not an arbitrary confidence score.

If visits and order value rise together, a component centred on [4 visits, $18] has an upward-tilted ellipse. That relationship is an off-diagonal covariance value.

covariance_type	What each component can learn	Practical consequence
full	Its own spread and rotation	Most flexible, but many parameters
tied	One shared full covariance	All components have the same shape
diag	Its own axis-aligned spread	More stable when rotation is unnecessary
spherical	Its own single radius	Round clusters; simplest option

A full covariance in d dimensions has d × (d + 1) / 2 distinct parameters. In 50 dimensions, that is 1,275 entries per component before means or weights. With only a few hundred rows, use diag, reduce dimension with PCA, or question whether a GMM is appropriate.

With full, diag, or tied covariance, converting dollars to cents is theoretically just an invertible rescaling: responsibilities and BIC rankings stay the same, although density units change. If they do not, inspect spherical covariance, fixed absolute reg_covar, convergence, and numerical conditioning.

Standardisation remains useful for conditioning and for making an absolute covariance floor meaningful. Fit its statistics on training data only. A log transform may better represent strongly right-skewed positive values.

EM: how the unknown groups become computable

The problem is circular:

• To estimate each Gaussian, we need to know which observations belong to it.
• To assign observations, we need the Gaussian parameters.

Expectation-Maximization, or EM, alternates between those calculations. It starts with provisional parameters and repeatedly improves them.

The E-step

Hold the means, covariances, and weights fixed. For every observation, calculate its responsibility for every component. For the 8-visit customer, the E-step produced approximately 0.293 for A and 0.707 for B.

The M-step

Hold the responsibilities fixed and re-estimate each component as if they were fractional counts.

Take four customers with weekly visits [3, 5, 9, 11]. Suppose the high-visit component H receives responsibilities:

[0.05, 0.20, 0.80, 0.95]

Its effective count is 2.00, so:

π_H = (0.05 + 0.20 + 0.80 + 0.95) / 4 = 0.50

Its new mean is:

μ_H = (0.05 × 3 + 0.20 × 5 + 0.80 × 9 + 0.95 × 11) / 2.00 = 9.4

The low component’s mean becomes 4.6. The weighted variance of each component is 4.24 visits squared, so both standard deviations become about 2.06.

In multiple dimensions, the weighted variance becomes a weighted covariance matrix. High-responsibility points pull strongly; low-responsibility points barely do.

EM repeats the E-step with the updated Gaussians, then another M-step. It stops when the log-likelihood or parameter updates barely change. Different initial parameters can lead to different local optima, so use multiple restarts.

The log-likelihood is the sum of the log mixture density for every training observation. Higher is better on the same data, though it does not guarantee a useful or generalisable segmentation.

Choosing the number of components

A GMM needs K, the number of components. A component is a statistical explanation, not necessarily a human category.

Fit several candidate values and compare an information criterion. BIC is:

BIC = -2 × log-likelihood + p × log(n)

Here, p is the number of fitted parameters and n is the number of observations. Lower BIC is better: the first term rewards fit and the second penalises complexity. AIC uses 2 × p instead, so it usually penalises complexity less.

Compare candidates on the same rows and preprocessing. Inspect stability across restarts and whether the groups support a useful decision. For predictive density, evaluate held-out log-likelihood. For segmentation, assess stability and actionability. A lower BIC does not make component labels semantically meaningful; it only indicates a better fit-complexity trade-off under the Gaussian assumptions.

A small working example

This code creates two customer-like populations, fits four candidate values of K, uses ten EM restarts, and selects the lowest BIC. It then finds the most ambiguous observation.

The exact BIC values depend on the installed numerical libraries.

import numpy as np

from sklearn.mixture import GaussianMixture

rng = np.random.default_rng(7)

X = np.vstack([
    rng.multivariate_normal(
        mean=[4, 18],
        cov=[[1.5, 0.8], [0.8, 4.0]],
        size=300,
    ),
    rng.multivariate_normal(
        mean=[8, 24],
        cov=[[4.0, -1.2], [-1.2, 2.5]],
        size=300,
    ),
])

models = {}
for k in range(1, 5):
    model = GaussianMixture(
        n_components=k,
        covariance_type="full",
        n_init=10,
        reg_covar=1e-6,
        random_state=0,
    )
    models[k] = model.fit(X)
    print(
        f"k={k}, BIC={model.bic(X):.1f}, "
        f"converged={model.converged_}"
    )

best_k = min(models, key=lambda k: models[k].bic(X))
gmm = models[best_k]

probabilities = gmm.predict_proba(X)
uncertain = np.argmin(np.max(probabilities, axis=1))

print("selected components:", best_k)
print(
    "most ambiguous probabilities:",
    np.round(probabilities[uncertain], 3),
)
print(
    "mixture density:",
    np.exp(gmm.score_samples(X[uncertain:uncertain + 1]))[0],
)

predict_proba returns one responsibility vector per row. score_samples returns the log of the overall mixture density. Density and membership differ: a point can be 90% assigned to one component while still being unusual under every component.

For anomaly scoring, use overall negative log density rather than the smallest membership probability. The latter only says that no component wins decisively.

GMM or something else?

Choose based on geometry and output:

Need or data shape	Usually start with	Why
Roughly round groups and hard labels	K-means	Fast squared-distance objective
Overlap probabilities and elliptical groups	GMM	Models density, uncertainty, and covariance
Curved groups or explicit noise points	DBSCAN	Follows density-connected shapes and marks noise
Nested groupings	Hierarchical clustering	A merge tree can be more useful than one partition

A GMM is a poor choice for crescents, rings, or long winding paths. Many components can approximate those shapes, but the resulting components may be arbitrary fragments. GMMs also cost more than k-means: full covariance has quadratic parameter growth in the number of features. In high dimensions, apparent confidence may be covariance overfitting. diag is less expressive but often a better baseline.

Failure modes you can see first

First symptom	Likely cause	Fix
Repeated fits give different groups or converged_ is false	Poor starts or different local optima	Scale features, increase n_init, and inspect several solutions
An ellipse collapses around one or two observations	Nearly singular covariance	Increase reg_covar, constrain covariance, reduce K, or use diag/tied
BIC keeps falling as K grows and tiny components appear	Extra components memorise tails or non-Gaussian shape	Check held-out density, component sizes, and another clustering family
Responsibilities are nearly all 0.000 or 1.000	Overconfidence, misspecification, or poor scaling	Check held-out likelihood, transformations, and covariance choice

The singular-covariance problem is fundamental. In an unconstrained mixture, a component can centre on one training point and shrink its covariance almost to zero. Its density, and therefore the likelihood, can become arbitrarily large. That is a mathematical loophole, not a microscopic customer segment.

A positive covariance floor such as reg_covar, explicit covariance constraints, or a model with a bounded likelihood prevents the collapse. More data and fewer components improve stability but do not, by themselves, bound the unconstrained mixture likelihood.

Quick check

Quick check

0/3

Q1A customer receives GMM responsibilities of 0.30 for component A and 0.70 for component B. What does that mean?

The customer belongs to two real groups at the same time, with a verified 30% and 70% chanceGiven the fitted model and this customer's features, component B contributes more posterior probability than component AThe customer is 70% closer to B under the k-means distanceThe model has discovered that B contains 70% of all customers

Check answer

Q2What happens during the M-step?

The model calculates each point's probability under the current componentsThe model chooses a new number of components without refittingThe model uses the current responsibilities as fractional counts to update each component's weight, mean, and covarianceThe model converts soft assignments into permanent labels

Check answer

Q3Transfer: your data forms two crescent-shaped groups, includes many noise points, and you do not need membership probabilities. Why might a GMM be a poor first choice?

A GMM can only process one featureA GMM assumes its components are Gaussian ellipses, so several components may fragment the crescents; a density method such as DBSCAN is a better first experimentBIC cannot be calculated for a GMMEM only works when every point has a known label

Check answer

Next

For broader clustering geometry, see K-means clustering. For visualising high-dimensional data, see t-SNE and UMAP.