Beyond Gaussian Mixtures

Every mixture model we've studied—Gaussian mixtures, Student-t mixtures, mixtures of experts, HMMs—requires specifying the number of components $K$ upfront. This presents a fundamental dilemma:

• Set $K$ too small: The model underfits, failing to capture distinct structure
• Set $K$ too large: The model overfits, splitting genuine clusters into artificial fragments
• Model selection: Running multiple models with different $K$ and using BIC/AIC is computationally expensive and feels inelegant

What if the model could learn the appropriate complexity from data? What if we could specify a prior that says "use as many components as the data warrants, but no more"?

Dirichlet Process Mixture Models (DPMMs) provide exactly this capability. By placing a Dirichlet Process prior over the space of distributions, DPMMs allow for a theoretically infinite number of mixture components, with the posterior naturally concentrating on a finite (and appropriate) number based on the data.

This "infinite mixture" might sound paradoxical—how can we compute with infinitely many components? The key insight is that, almost surely, only finitely many components are ever used, and we can perform inference without explicitly representing the unused components.

Before introducing the Dirichlet Process, we need to thoroughly understand the Dirichlet distribution, which generalizes the Beta distribution to multiple dimensions.

Definition

The Dirichlet distribution $\text{Dir}(\boldsymbol{\alpha})$ with concentration parameters $\boldsymbol{\alpha} = (\alpha_1, \ldots, \alpha_K)$ is a distribution over the $(K-1)$-simplex—the set of vectors $\boldsymbol{\pi} = (\pi_1, \ldots, \pi_K)$ satisfying $\pi_k \geq 0$ and $\sum_k \pi_k = 1$.

Density: $$p(\boldsymbol{\pi} | \boldsymbol{\alpha}) = \frac{\Gamma(\sum_{k=1}^K \alpha_k)}{\prod_{k=1}^K \Gamma(\alpha_k)} \prod_{k=1}^K \pi_k^{\alpha_k - 1}$$

Key properties:

• $\mathbb{E}[\pi_k] = \alpha_k / \alpha_0$ where $\alpha_0 = \sum_k \alpha_k$
• $\text{Var}(\pi_k) = \frac{\alpha_k(\alpha_0 - \alpha_k)}{\alpha_0^2(\alpha_0 + 1)}$
• As $\alpha_0 \to \infty$, variance shrinks to zero

Interpretation of Concentration Parameters

The parameters $\alpha_k$ control the expected proportion and variability:

Symmetric case $\alpha_k = \alpha$ for all $k$:

• $\alpha = 1$: Uniform over the simplex
• $\alpha < 1$: Mass concentrates near vertices (sparse: most weight on few components)
• $\alpha > 1$: Mass concentrates near center (uniform: weight spread across components)

Conjugacy with Multinomial

The Dirichlet is conjugate to the Categorical/Multinomial:

$$\boldsymbol{\pi} \sim \text{Dir}(\boldsymbol{\alpha})$$ $$x_1, \ldots, x_N | \boldsymbol{\pi} \sim \text{Categorical}(\boldsymbol{\pi})$$ $$\boldsymbol{\pi} | x_1, \ldots, x_N \sim \text{Dir}(\boldsymbol{\alpha} + \mathbf{n})$$

where $n_k = \sum_i \mathbb{1}[x_i = k]$ counts observations in category $k$.

This closed-form update is why the Dirichlet is foundational for Bayesian mixture models.

The Dirichlet Process (DP) is a distribution over probability distributions—a "distribution over distributions." It extends the Dirichlet distribution from finite to infinite dimensions.

Formal Definition

A random probability measure $G$ follows a Dirichlet Process with base distribution $H$ and concentration parameter $\alpha$, written $G \sim \text{DP}(\alpha, H)$, if for any finite measurable partition $(A_1, \ldots, A_K)$ of the sample space:

$$(G(A_1), \ldots, G(A_K)) \sim \text{Dir}(\alpha H(A_1), \ldots, \alpha H(A_K))$$

Interpretation:

• $H$ is the base distribution—the expected value $\mathbb{E}[G] = H$
• $\alpha$ is the concentration parameter—controls how close $G$ is to $H$
• As $\alpha \to \infty$, $G \to H$ (the prior dominates)
• As $\alpha \to 0$, $G$ concentrates on few atoms

Key Property: Discreteness

A crucial property of the DP is that draws $G \sim \text{DP}(\alpha, H)$ are almost surely discrete, even when $H$ is continuous. This means $G$ can be written as:

$$G = \sum_{k=1}^\infty \pi_k \delta_{\theta_k}$$

where:

• $\delta_{\theta}$ is a point mass at $\theta$
• $\theta_k \sim H$ are the locations (drawn from base distribution)
• $\pi_k$ are the weights (summing to 1)

This decomposition is made explicit by the stick-breaking construction.

Why This Matters for Mixtures

In a standard mixture model: $$p(x) = \sum_{k=1}^K \pi_k f(x | \theta_k)$$

With a DP prior on the mixing distribution: $$G \sim \text{DP}(\alpha, H)$$ $$\theta_i | G \sim G$$ $$x_i | \theta_i \sim f(x | \theta_i)$$

Since $G$ is discrete, multiple observations can share the same $\theta$ value—they belong to the same cluster. But the number of distinct $\theta$ values (clusters) is random and grows with data.

The stick-breaking construction (Sethuraman, 1994) provides an explicit, constructive representation of draws from the Dirichlet Process.

The Procedure

To sample $G \sim \text{DP}(\alpha, H)$:

1. Generate stick-breaking weights: $$\beta_k \sim \text{Beta}(1, \alpha), \quad k = 1, 2, \ldots$$ $$\pi_k = \beta_k \prod_{j=1}^{k-1} (1 - \beta_j)$$

2. Generate atoms: $$\theta_k \sim H, \quad k = 1, 2, \ldots$$

3. The random measure: $$G = \sum_{k=1}^\infty \pi_k \delta_{\theta_k}$$

Stick-Breaking Intuition

Imagine a stick of unit length representing total probability mass:

1. Break off proportion $\beta_1$ from the remaining stick → weight $\pi_1 = \beta_1$
2. Break off proportion $\beta_2$ from what remains → weight $\pi_2 = \beta_2(1-\beta_1)$
3. Continue indefinitely...

Since each $\beta_k < 1$ (almost surely), the remaining stick shrinks exponentially, ensuring weights sum to 1 despite the infinite process.

Properties of Stick-Breaking Weights

Expected weight: $$\mathbb{E}[\pi_k] = \frac{1}{\alpha + 1} \left(\frac{\alpha}{\alpha + 1}\right)^{k-1}$$

Weights decay geometrically with rate $\alpha/(\alpha+1)$.

Effect of $\alpha$:

• Small $\alpha$ (e.g., 0.1): $\beta_k \approx 1$ often, so $\pi_1$ captures most mass. Few clusters.
• Large $\alpha$ (e.g., 100): $\beta_k \approx 0$ often, so mass spreads across many components.

Truncation for Computation

For practical inference, truncate to $K$ components: $$G_K = \sum_{k=1}^{K-1} \pi_k \delta_{\theta_k} + \left(1 - \sum_{k=1}^{K-1} \pi_k\right) \delta_{\theta_K}$$

The last component absorbs all remaining mass. With $K$ sufficiently large, this approximation is excellent: $$\mathbb{E}\left[\sum_{k=K+1}^\infty \pi_k\right] = \left(\frac{\alpha}{\alpha + 1}\right)^K$$

For $\alpha = 1$ and $K = 50$, the truncation error is $\approx 10^{-15}$.

The Chinese Restaurant Process (CRP) provides an alternative, sequential characterization of the Dirichlet Process that's particularly useful for understanding clustering behavior and deriving inference algorithms.

The Metaphor

Imagine a restaurant with infinitely many tables (clusters), each with unlimited capacity. Customers (data points) arrive sequentially and choose where to sit:

Customer 1: Sits at the first table.

Customer $n+1$: Given $n$ customers already seated:

• Join existing table $k$ (with $n_k$ customers) with probability $\frac{n_k}{n + \alpha}$
• Start a new table with probability $\frac{\alpha}{n + \alpha}$

Mathematical Form

Let $z_i$ be the table assignment for customer $i$. The conditional distribution is:

$$p(z_{n+1} = k | z_1, \ldots, z_n) = \begin{cases} \frac{n_k}{n + \alpha} & \text{if } k \text{ is an existing table} \ \frac{\alpha}{n + \alpha} & \text{if } k \text{ is a new table} \end{cases}$$

where $n_k = \sum_{i=1}^n \mathbb{1}[z_i = k]$ is the number of customers at table $k$.

Key Properties

1. Rich-get-richer: Popular tables attract more customers (preferential attachment). This leads to power-law table size distributions.

2. New tables always possible: Non-zero probability of a new table ensures cluster count grows with data.

3. Expected number of tables: For $n$ customers: $$\mathbb{E}[K_n] \approx \alpha \log(n/\alpha + 1)$$

Logarithmic growth—much slower than linear.

4. Exchangeability: The joint distribution $p(z_1, \ldots, z_n)$ is exchangeable—invariant to customer ordering. The sequential construction is for conceptual convenience.

Pólya Urn Representation

An equivalent view uses urns and balls:

1. Start with urn containing $\alpha$ "special" balls
2. To seat customer $n+1$:
   • Draw a ball uniformly at random
   • If special ball: new table, add regular ball for that table
   • If regular ball of color $k$: join table $k$, return ball, add another of same color

The DPMM combines the Dirichlet Process prior with a mixture model likelihood, enabling nonparametric clustering.

Generative Model

Stick-breaking formulation: $$\beta_k \sim \text{Beta}(1, \alpha), \quad \pi_k = \beta_k \prod_{j < k}(1 - \beta_j)$$ $$\theta_k \sim H \quad \text{(component parameters)}$$ $$z_i \sim \text{Categorical}(\boldsymbol{\pi}) \quad \text{(cluster assignment)}$$ $$x_i \sim F(\theta_{z_i}) \quad \text{(observation)}$$

CRP formulation: $$z_i | z_1, \ldots, z_{i-1} \sim \text{CRP}(\alpha)$$ $$\theta_k \sim H \quad \text{for each unique cluster } k$$ $$x_i | z_i, {\theta_k} \sim F(\theta_{z_i})$$

Gaussian DPMM Example

With Gaussian observations and Normal-Inverse-Wishart base: $$H = \text{NIW}(\mu_0, \kappa_0, u_0, \boldsymbol{\Psi}_0)$$ $$\theta_k = (\boldsymbol{\mu}_k, \boldsymbol{\Sigma}_k) \sim H$$ $$x_i | z_i = k \sim \mathcal{N}(\boldsymbol{\mu}_k, \boldsymbol{\Sigma}_k)$$

The NIW is conjugate to the Gaussian likelihood, enabling closed-form updates.

Predictive Distribution

A beautiful property of DPMMs: after observing data, the predictive distribution for a new point has a simple form.

Using the CRP, the predictive assignment for $x_{n+1}$:

$$p(z_{n+1} | z_{1:n}, x_{1:n+1}) \propto \begin{cases} n_k \cdot p(x_{n+1} | x_{\text{in cluster } k}) & \text{existing cluster } k \ \alpha \cdot p(x_{n+1} | H) & \text{new cluster} \end{cases}$$

For conjugate models, $p(x_{n+1} | x_{\text{in cluster } k})$ integrates out $\theta_k$ analytically (the posterior predictive).

Comparison to Finite Mixtures

Posterior inference in DPMMs is intractable analytically, requiring sampling or variational approximations.

Collapsed Gibbs Sampling

The most common approach marginalizes out component parameters $\theta$ and samples cluster assignments $z$ directly.

Algorithm (Neal's Algorithm 2):

For each data point $i$:

1. Remove $x_i$ from its current cluster, decrementing $n_{z_i}$
2. If cluster becomes empty, remove it
3. Sample new assignment from: $$p(z_i = k | z_{-i}, \mathbf{x}, \alpha) \propto \begin{cases} n_{-i,k} \cdot p(x_i | \mathbf{x}_{-i, k}) & k \text{ existing} \ \alpha \cdot p(x_i | H) & k = \text{new} \end{cases}$$
4. Add $x_i$ to sampled cluster

Here $z_{-i}$ denotes assignments excluding $i$, and $\mathbf{x}_{-i,k}$ are observations in cluster $k$ excluding $x_i$.

Key: $p(x_i | \mathbf{x}_{-i,k})$ is the posterior predictive, which is analytical for conjugate models.

Variational Inference

For large datasets, MCMC can be slow. Variational inference provides a faster alternative.

Truncated stick-breaking VI:

1. Truncate to $K$ components (large enough)
2. Approximate posterior: $q(\boldsymbol{\beta}, \boldsymbol{\theta}, \mathbf{z}) = \prod_k q(\beta_k) q(\theta_k) \prod_i q(z_i)$
3. Optimize ELBO via coordinate ascent

Stochastic Variational Inference (SVI):

• Process mini-batches of data
• Scales to millions of points
• Implemented in libraries like scikit-learn and PyMC

Applying DPMMs effectively requires careful attention to several practical issues.

Choosing the Concentration Parameter α

Option 1: Fix based on expected clusters

• Use formula $\mathbb{E}[K] \approx \alpha \log(N/\alpha + 1)$
• If you expect ~10 clusters for 1000 points, solve for α

Option 2: Hierarchical prior

• Place $\alpha \sim \text{Gamma}(a, b)$ and sample it
• Common choice: $a = 1$, $b = 1$ (vague prior)

Option 3: Cross-validation

• Evaluate held-out likelihood for different α values
• Select α that maximizes predictive performance

Base Distribution Selection

The base distribution $H$ should:

• Cover the support of possible cluster parameters
• Be conjugate to the likelihood for computational efficiency
• Encode reasonable prior beliefs

Common choices:

• Gaussian likelihood → NIW base
• Categorical likelihood → Dirichlet base
• Poisson likelihood → Gamma base

This page has provided a comprehensive treatment of Dirichlet Process Mixture Models, the foundational nonparametric Bayesian approach to clustering.

What's Next: Semi-Parametric Methods

DPMMs represent fully nonparametric density estimation—flexible but computationally demanding. The final page of this module explores semi-parametric methods that strike a balance: using parametric models locally while allowing nonparametric flexibility globally. This includes mixture density networks, normalizing flows for mixture models, and copula-based approaches.