Gaussian Mixture Models

In statistical modeling and machine learning, we frequently encounter data that cannot be adequately described by a single parametric distribution. Consider a dataset of human heights that includes both adult males and adult females—a single Gaussian distribution would fail to capture the bimodal nature of the data. Or imagine customer purchase amounts in a retail dataset, where different customer segments exhibit fundamentally different spending behaviors.

Mixture models provide an elegant and principled solution to this challenge. Rather than forcing all data through a single distributional lens, mixture models posit that the observed data arises from multiple underlying subpopulations, each characterized by its own parametric distribution. The overall data distribution is then a weighted combination—a mixture—of these component distributions.

This page develops the complete mathematical formulation of Gaussian Mixture Models (GMMs), the most widely used instance of mixture modeling. We will derive every component of the model from first principles, ensuring you understand not just what the equations say, but why they take the form they do.

Before constructing mixture models, we must establish a rigorous understanding of the Gaussian distribution itself. The univariate Gaussian (or Normal) distribution is parameterized by two quantities: the mean $\mu$ (location parameter) and the variance $\sigma^2$ (scale parameter). Its probability density function (pdf) is:

Dissecting the components:

The normalization constant $\frac{1}{\sqrt{2\pi\sigma^2}}$ ensures that the pdf integrates to 1 over the entire real line. This is a requirement for any valid probability density function. Notice that larger variance $\sigma^2$ leads to a smaller normalization constant—because the probability mass is spread over a wider range, the peak height must decrease.

The exponential term $\exp\left(-\frac{(x - \mu)^2}{2\sigma^2}\right)$ captures the shape of the distribution. The squared term $(x - \mu)^2$ measures the distance from the mean, ensuring symmetry around $\mu$. The division by $2\sigma^2$ scales this distance relative to the spread—points many standard deviations from the mean receive exponentially lower density.

The maximum density occurs at $x = \mu$, where the exponent is zero and the exponential term equals 1. As $|x - \mu|$ increases, the density decreases rapidly due to the negative quadratic in the exponent.

Real-world data typically involves multiple dimensions. Customer features might include age, income, and purchase frequency. Image data involves thousands of pixel intensities. To handle such data, we extend the Gaussian to multiple dimensions.

For a $D$-dimensional random vector $\mathbf{x} = (x_1, x_2, \ldots, x_D)^\top$, the multivariate Gaussian distribution is parameterized by:

• Mean vector $\boldsymbol{\mu} \in \mathbb{R}^D$: The center of the distribution in $D$-dimensional space
• Covariance matrix $\boldsymbol{\Sigma} \in \mathbb{R}^{D \times D}$: A symmetric, positive definite matrix capturing variance and pairwise correlations

Understanding each component:

The normalization constant $(2\pi)^{D/2}|\boldsymbol{\Sigma}|^{1/2}$ generalizes the univariate case. The $D/2$ power accounts for $D$ independent dimensions, while the determinant $|\boldsymbol{\Sigma}|$ measures the "volume" of the distribution's spread. Larger determinant means more spread, requiring a smaller peak height.

The Mahalanobis distance $(\mathbf{x} - \boldsymbol{\mu})^\top \boldsymbol{\Sigma}^{-1} (\mathbf{x} - \boldsymbol{\mu})$ is a generalized squared distance that accounts for:

• Different variances along each dimension
• Correlations between dimensions

This quadratic form defines elliptical contours of constant density. Points on the same ellipse have equal probability density.

The precision matrix $\boldsymbol{\Sigma}^{-1}$, inverse of the covariance, appears rather than $\boldsymbol{\Sigma}$ itself. This is computationally important—we often work with precision directly in optimization.

The parameter count matters for mixture models. A full covariance matrix in $D$ dimensions has $D(D+1)/2$ free parameters due to symmetry. For a mixture with $K$ components, this means $K \cdot D(D+1)/2$ covariance parameters, which can quickly become prohibitive. This motivates restricted forms:

Spherical covariance (isotropic): All dimensions have equal variance, no correlation. This produces spherical clusters.

Diagonal covariance: Each dimension has its own variance, but no correlation between dimensions. This produces axis-aligned elliptical clusters.

Full covariance: Arbitrary orientation and scaling. Most flexible but most parameters.

Tied covariance: All components share the same covariance matrix, differing only in means. Reduces parameters significantly.

Now we construct the Gaussian Mixture Model by combining multiple Gaussian components. The key insight is that complex, multimodal distributions can often be approximated as weighted sums of simpler unimodal distributions.

A Gaussian Mixture Model with $K$ components is defined by:

1. $K$ component densities: Each component $k \in {1, 2, \ldots, K}$ is a multivariate Gaussian with its own parameters $(\boldsymbol{\mu}_k, \boldsymbol{\Sigma}_k)$

2. $K$ mixing coefficients: Non-negative weights $\pi_1, \pi_2, \ldots, \pi_K$ that sum to one: $\sum_{k=1}^{K} \pi_k = 1$

Why this is a valid probability density:

For $p(\mathbf{x})$ to be a valid pdf, it must:

1. Be non-negative everywhere: $p(\mathbf{x}) \geq 0$ for all $\mathbf{x}$
2. Integrate to one: $\int p(\mathbf{x}) d\mathbf{x} = 1$

Non-negativity is immediate: each $\pi_k \geq 0$ and each $\mathcal{N}(\mathbf{x} \mid \boldsymbol{\mu}_k, \boldsymbol{\Sigma}_k) \geq 0$, so their product and sum are non-negative.

Integration to one follows from: $$\int p(\mathbf{x}) d\mathbf{x} = \int \sum_{k=1}^{K} \pi_k , \mathcal{N}(\mathbf{x} \mid \boldsymbol{\mu}_k, \boldsymbol{\Sigma}k) d\mathbf{x}$$ $$= \sum{k=1}^{K} \pi_k \int \mathcal{N}(\mathbf{x} \mid \boldsymbol{\mu}_k, \boldsymbol{\Sigma}k) d\mathbf{x} = \sum{k=1}^{K} \pi_k \cdot 1 = 1$$

The integral and sum exchange by linearity, and each Gaussian integrates to 1.

The complete parameter set:

Collecting all parameters, a GMM with $K$ components in $D$ dimensions has:

$$\boldsymbol{\theta} = {\pi_1, \ldots, \pi_K, \boldsymbol{\mu}_1, \ldots, \boldsymbol{\mu}_K, \boldsymbol{\Sigma}_1, \ldots, \boldsymbol{\Sigma}_K}$$

Parameter count for full covariance:

• Mixing coefficients: $K - 1$ (due to sum-to-one constraint)
• Means: $K \cdot D$
• Covariances: $K \cdot D(D+1)/2$

Total: $K - 1 + KD + K \cdot D(D+1)/2 = K(1 + D + D(D+1)/2) - 1$ parameters

For example, with $K = 5$ components in $D = 10$ dimensions with full covariance:

• $5 - 1 = 4$ mixing coefficient parameters
• $5 \times 10 = 50$ mean parameters
• $5 \times 55 = 275$ covariance parameters
• Total: 329 parameters

This illustrates why simpler covariance structures are often necessary in practice.

One of the most powerful ways to understand probabilistic models is through their generative story—a step-by-step description of how the observed data could have been produced.

For a Gaussian Mixture Model, the generative process for a single observation is:

Step 1: Sample the component assignment Draw a component index $z$ from a Categorical distribution with probabilities $(\pi_1, \ldots, \pi_K)$: $$z \sim \text{Categorical}(\pi_1, \ldots, \pi_K)$$

Step 2: Sample the observation from the chosen component Given $z = k$, draw the observation $\mathbf{x}$ from the $k$-th Gaussian: $$\mathbf{x} \mid z = k \sim \mathcal{N}(\boldsymbol{\mu}_k, \boldsymbol{\Sigma}_k)$$

To generate a dataset of $N$ observations:

Repeat the two-step process $N$ times independently. Each observation $\mathbf{x}_n$ is associated with an unobserved (latent) component assignment $z_n$. The complete generative process produces:

• Observations: ${\mathbf{x}_1, \mathbf{x}_2, \ldots, \mathbf{x}_N}$ (what we observe)
• Assignments: ${z_1, z_2, \ldots, z_N}$ (hidden/latent)

The fundamental challenge:

In practice, we observe only the $\mathbf{x}_n$ values—the component assignments $z_n$ are latent (hidden). We cannot directly count which observations belong to which component. This is precisely what makes mixture model inference challenging and leads to the Expectation-Maximization algorithm (covered in Module 4).

We now show that the generative story rigorously implies the mixture density formula. This derivation uses fundamental probability rules and deepens understanding of the model structure.

Starting point: Joint distribution of $\mathbf{x}$ and $z$

The generative story specifies:

• Prior on component: $P(z = k) = \pi_k$
• Likelihood given component: $p(\mathbf{x} \mid z = k) = \mathcal{N}(\mathbf{x} \mid \boldsymbol{\mu}_k, \boldsymbol{\Sigma}_k)$

By the chain rule of probability: $$p(\mathbf{x}, z = k) = P(z = k) \cdot p(\mathbf{x} \mid z = k) = \pi_k , \mathcal{N}(\mathbf{x} \mid \boldsymbol{\mu}_k, \boldsymbol{\Sigma}_k)$$

The law of total probability in action:

The marginalization step is an application of the law of total probability: $$p(\mathbf{x}) = \sum_{k} P(z = k) \cdot p(\mathbf{x} \mid z = k)$$

This says: the probability of observing $\mathbf{x}$ is the sum over all ways $\mathbf{x}$ could have been generated—weighting each component's contribution by the probability of selecting that component.

Posterior responsibility:

Conversely, given an observation $\mathbf{x}$, we can compute the posterior probability that it came from component $k$ using Bayes' theorem: $$P(z = k \mid \mathbf{x}) = \frac{P(z = k) \cdot p(\mathbf{x} \mid z = k)}{p(\mathbf{x})} = \frac{\pi_k , \mathcal{N}(\mathbf{x} \mid \boldsymbol{\mu}_k, \boldsymbol{\Sigma}k)}{\sum{j=1}^{K} \pi_j , \mathcal{N}(\mathbf{x} \mid \boldsymbol{\mu}_j, \boldsymbol{\Sigma}_j)}$$

This quantity, often denoted $\gamma_{nk}$ for observation $n$ and component $k$, is called the responsibility—how responsible component $k$ is for generating observation $\mathbf{x}_n$. It plays a central role in the EM algorithm.

Given observed data $\mathbf{X} = {\mathbf{x}_1, \ldots, \mathbf{x}_N}$, we want to estimate the GMM parameters $\boldsymbol{\theta} = {\pi_k, \boldsymbol{\mu}_k, \boldsymbol{\Sigma}k}{k=1}^{K}$. The maximum likelihood approach seeks parameters that maximize the probability of observing the data.

Assuming independence:

If observations are independent and identically distributed (i.i.d.) given the parameters, the joint likelihood is: $$L(\boldsymbol{\theta}) = p(\mathbf{X} \mid \boldsymbol{\theta}) = \prod_{n=1}^{N} p(\mathbf{x}n \mid \boldsymbol{\theta}) = \prod{n=1}^{N} \sum_{k=1}^{K} \pi_k , \mathcal{N}(\mathbf{x}_n \mid \boldsymbol{\mu}_k, \boldsymbol{\Sigma}_k)$$

Why the GMM log-likelihood is difficult to optimize:

Note the structure of the log-likelihood. The logarithm is outside the summation over components. This is the source of all optimization difficulty:

$$\ell(\boldsymbol{\theta}) = \sum_{n=1}^{N} \log \left( \sum_{k=1}^{K} \pi_k , \mathcal{N}(\mathbf{x}_n \mid \boldsymbol{\mu}_k, \boldsymbol{\Sigma}_k) \right)$$

Comparison with single Gaussian:

For a single Gaussian, the log-likelihood is: $$\ell(\boldsymbol{\mu}, \boldsymbol{\Sigma}) = \sum_{n=1}^{N} \log \mathcal{N}(\mathbf{x}n \mid \boldsymbol{\mu}, \boldsymbol{\Sigma})$$ $$= -\frac{N}{2} \log|\boldsymbol{\Sigma}| - \frac{1}{2} \sum{n=1}^{N} (\mathbf{x}_n - \boldsymbol{\mu})^\top \boldsymbol{\Sigma}^{-1} (\mathbf{x}_n - \boldsymbol{\mu}) + \text{const}$$

This is a concave function of the parameters (after suitable reparameterization), and the maximum has a closed-form solution: $\hat{\boldsymbol{\mu}} = \frac{1}{N}\sum_n \mathbf{x}_n$ and $\hat{\boldsymbol{\Sigma}} = \frac{1}{N}\sum_n (\mathbf{x}_n - \hat{\boldsymbol{\mu}})(\mathbf{x}_n - \hat{\boldsymbol{\mu}})^\top$.

The mixture complication:

For the GMM, the log-of-sum structure prevents closed-form solutions. Setting derivatives to zero yields equations where each parameter appears on both sides, entangled with all other parameters. No algebraic rearrangement isolates them.

This is why we need iterative algorithms like Expectation-Maximization (EM)—the subject of Module 4.

Understanding mixture densities geometrically builds intuition for how component parameters affect the overall distribution. Let's examine how different configurations produce different shapes.

One-dimensional mixtures:

In 1D, mixture densities appear as weighted sums of bell curves. The resulting shape exhibits:

• Multiple modes when components are well-separated
• Shoulders or asymmetry when components overlap partially
• Approximate unimodality when components overlap heavily

Two-dimensional mixtures:

In 2D, we can visualize with contour plots or 3D surface plots. Each Gaussian component produces elliptical contours (circles for spherical covariance). The mixture density shows:

• Complex, non-convex level sets
• Ridges connecting nearby modes
• Valleys between well-separated components

Key observations from visualization:

1. Superposition principle: The mixture density at any point is the sum of individual Gaussian densities. Where components overlap, densities add, creating ridges or elevated regions between modes.

2. Non-elliptical contours: Unlike a single Gaussian, mixture contours can have complex, non-convex shapes. Contours at lower density levels may encircle multiple modes.

3. Weight effects: Components with larger $\pi_k$ contribute more to the overall density. A component with small weight barely affects the density shape, even if its Gaussian is well-separated.

4. Covariance effects: The orientation and eccentricity of each component's ellipse affects how components "blend" together. Aligned components merge smoothly; perpendicular ones create more complex intersections.

5. Separation effects: Well-separated components produce clearly multimodal distributions. Overlapping components create ambiguity in component assignment.

We have developed the complete mathematical formulation of Gaussian Mixture Models from first principles. Let's consolidate the key concepts:

What comes next:

This formulation page establishes the probabilistic framework. The subsequent pages in this module build upon this foundation:

• Page 1: Latent Variables — Deepens the understanding of $z$ as a latent variable, exploring the complete-data likelihood and why latent variables make GMMs both powerful and challenging.

• Page 2: Marginal Likelihood — Develops the marginal likelihood perspective for model comparison and its connection to Bayesian inference.

• Page 3: Identifiability — Addresses when GMM parameters are uniquely determined and the label-switching symmetry.

• Page 4: Model Selection — Covers methods for determining the optimal number of components $K$.

The EM algorithm for fitting GMMs is covered in detail in Module 4: Expectation-Maximization Algorithm.