Bayesian Nonparametrics

Overview & Motivation

In parametric statistics, we choose a model family — say, Gaussian distributions $\mathcal{N}(\mu, \sigma^2)$ — and reduce inference to estimating a fixed, finite-dimensional parameter $\theta = (\mu, \sigma^2) \in \mathbb{R}^2$. This works beautifully when the model is well-specified. But what if the true data-generating process is multimodal, heavy-tailed, or otherwise poorly captured by any finite-dimensional family?

The PAC learning framework gave us one answer: control model complexity through the VC dimension or Rademacher complexity, using structural risk minimization (SRM) to balance approximation error and estimation error. Bayesian nonparametrics offers a fundamentally different approach: place a prior directly on an infinite-dimensional parameter space and let the effective complexity of the posterior grow with the data.

The distinction between the two paradigms is worth making precise:

• Parametric models: fix the number of parameters a priori (e.g., fit a mixture of $K = 3$ Gaussians).
• Nonparametric models: let the effective number of parameters grow with $n$ (e.g., the number of mixture components adapts to the data).

This naming is somewhat misleading — “nonparametric” models have more parameters than parametric ones, not fewer. A better name might be “infinite-parametric,” but the convention is firmly established.

We’ll develop three canonical nonparametric models in this topic:

1. The Dirichlet Process — a prior on probability measures, used for clustering and density estimation.
2. The Gaussian Process — a prior on functions, used for regression and classification.
3. The Indian Buffet Process — a prior on binary matrices, used for latent feature models.

Each places a prior on an infinite-dimensional object (a measure, a function, a binary matrix with infinitely many columns), yet admits tractable posterior inference through clever constructive representations.

---

The Dirichlet Distribution

The Dirichlet process is the infinite-dimensional generalization of the Dirichlet distribution, so we begin by reviewing the finite-dimensional case carefully.

The concentration parameter $\alpha_0 = \sum_k \alpha_k$ controls how concentrated the distribution is around the base measure $\mathbf{b} = \boldsymbol{\alpha}/\alpha_0$:

• When $\alpha_0 \gg 1$, draws are concentrated near $\mathbf{b}$ (low variance).
• When $\alpha_0 \ll 1$, draws are concentrated near the vertices of the simplex (sparse, winner-take-all).
• When $\alpha_0 = 1$ and $\alpha_k = 1/K$ for all $k$, draws are approximately uniform on the simplex.

---

The Dirichlet Process

The Dirichlet process, introduced by Ferguson (1973), is the cornerstone of Bayesian nonparametrics. It is a distribution over probability distributions — a “prior over priors” — that generalizes the Dirichlet distribution to infinite-dimensional spaces.

This definition is elegant but requires verification that such an object exists — the condition must be self-consistent across all possible partitions.

The two parameters of the DP have clear roles:

• Base measure $G_0$: the “prior guess” at what $G$ looks like. $\mathbb{E}[G(A)] = G_0(A)$ for every measurable $A$ — draws from the DP are centered around $G_0$.
• Concentration parameter $\alpha$: controls how close $G$ is to $G_0$. As $\alpha \to \infty$, $G \to G_0$ in distribution. As $\alpha \to 0$, $G$ concentrates on a single atom.

The most surprising — and practically important — property of the Dirichlet process is its almost sure discreteness.

---

Constructive Representations

Ferguson’s definition (Definition 2) is clean but non-constructive — it tells us the finite-dimensional marginals without directly telling us how to sample from the DP. Three equivalent constructive representations fill this gap, each offering different computational and conceptual advantages.

The Stick-Breaking Construction

The name is vivid: imagine a stick of length 1. Break off a fraction $V_1$ (the first weight $w_1 = V_1$). From the remaining piece of length $1 - V_1$, break off a fraction $V_2$ (giving $w_2 = V_2(1 - V_1)$). Continue ad infinitum. The process almost surely exhausts the stick: $\sum_k w_k = 1$ with probability one.

In practice, we truncate the stick-breaking construction at $K$ components, which gives an excellent approximation when $K$ is large enough:

def stick_breaking_sample(alpha, G0_sampler, K=200):
    """Sample a truncated DP via stick-breaking."""
    V = rng.beta(1, alpha, size=K)
    w = np.zeros(K)
    w[0] = V[0]
    for k in range(1, K):
        w[k] = V[k] * np.prod(1 - V[:k])
    atoms = G0_sampler(K)
    return w, atoms

G0_sampler = lambda K: rng.normal(0, 1, K)

Stick-Breaking Construction Explorer

α (concentration): 1.0

K (truncation): 20

Resample

With α = 1.0, the DP draw is a discrete distribution over K = 20 atoms. Remaining stick mass: 0.0%.

The Chinese Restaurant Process

The CRP provides an intuitive simulation algorithm:

alpha_crp = 2.0
n_customers = 50
tables = []       # table sizes
assignments = []  # which table each customer sits at

for i in range(n_customers):
    if len(tables) == 0:
        tables.append(1)
        assignments.append(0)
    else:
        probs = np.array(tables + [alpha_crp]) / (i + alpha_crp)
        choice = rng.choice(len(probs), p=probs)
        if choice == len(tables):   # new table
            tables.append(1)
            assignments.append(len(tables) - 1)
        else:                        # existing table
            tables[choice] += 1
            assignments.append(choice)

The Pólya Urn Scheme

---

Posterior Inference

One of the most elegant properties of the Dirichlet process is its conjugacy: the posterior of a DP prior given i.i.d. observations is again a DP, with parameters updated in a natural way.

DP Posterior Updating

Click on the plot to add data points (max 20)

α = 1.0n = 0 observationsClear

---

Dirichlet Process Mixture Models

The DP by itself generates discrete distributions, but real data is often continuous. The Dirichlet process mixture model (DPMM) solves this by using the DP as a mixing distribution: each observation is drawn from a kernel (e.g., Gaussian) centered at a DP-sampled atom. This produces a countable mixture with an unknown number of components.

The generative process via stick-breaking:

1. Draw weights: $V_k \sim \text{Beta}(1, \alpha)$, set $w_k = V_k \prod_{j<k}(1-V_j)$.
2. Draw atoms: $\theta_k \sim G_0$.
3. For each observation $i$: assign to component $z_i$ with $\Pr[z_i = k] = w_k$, then draw $x_i \sim F(\cdot \mid \theta_{z_i})$.

Gibbs sampling for DPMMs. The collapsed Gibbs sampler (Neal, Algorithm 3) iterates over observations, resampling each cluster assignment $z_i$ from its full conditional:

$\Pr[z_i = k \mid z_{-i}, x_1, \ldots, x_n] \propto \begin{cases} n_{k,-i} \cdot F(x_i \mid \theta_k) & \text{existing cluster } k, \\ \alpha \cdot \int F(x_i \mid \theta) \, dG_0(\theta) & \text{new cluster.} \end{cases}$

When the kernel $F$ and base measure $G_0$ are conjugate (e.g., Gaussian-NIW), the marginal likelihood $\int F(x_i \mid \theta) dG_0(\theta)$ and the updated parameter $\theta_k$ given all observations assigned to cluster $k$ have closed-form expressions.

def run_dpmm_gibbs(X, alpha, n_iter=50, sigma2=0.5, mu0=0, sigma02=10):
    """Collapsed Gibbs sampler for a univariate Gaussian DPMM."""
    n = len(X)
    z = rng.integers(0, 3, size=n)  # random initial assignments
    for iteration in range(n_iter):
        for i in range(n):
            # Remove observation i from its cluster
            clusters = {}
            for j in range(n):
                if j == i: continue
                clusters.setdefault(z[j], []).append(j)
            probs, cluster_ids = [], []
            for k, members in clusters.items():
                n_k = len(members)
                x_bar = np.mean(X[members])
                precision_post = 1/sigma02 + n_k/sigma2
                mu_post = (mu0/sigma02 + n_k*x_bar/sigma2) / precision_post
                sigma2_pred = sigma2 + 1/precision_post
                log_prob = np.log(n_k) + norm.logpdf(X[i], mu_post, np.sqrt(sigma2_pred))
                probs.append(np.exp(log_prob))
                cluster_ids.append(k)
            # New cluster probability
            log_new = np.log(alpha) + norm.logpdf(X[i], mu0, np.sqrt(sigma2 + sigma02))
            probs.append(np.exp(log_new))
            cluster_ids.append(max(cluster_ids, default=-1) + 1)
            probs = np.array(probs) / sum(probs)
            z[i] = cluster_ids[rng.choice(len(probs), p=probs)]
    return z

---

Gaussian Processes

While the Dirichlet process provides a nonparametric prior on probability measures, the Gaussian process provides a nonparametric prior on functions. In the Bayesian nonparametric view, a GP is a prior on an infinite-dimensional function space that admits tractable finite-dimensional marginals.

The real power of GPs lies in the closed-form posterior — conditioning on observed data is just Gaussian conditioning.

In practice, the matrix inverse $(\mathbf{K} + \sigma_n^2 \mathbf{I})^{-1}$ is computed via the Cholesky decomposition for numerical stability:

def rbf_kernel(X1, X2, length_scale=1.0, signal_var=1.0):
    """Squared exponential (RBF) kernel."""
    sqdist = np.sum(X1**2, 1).reshape(-1, 1) + np.sum(X2**2, 1) - 2 * X1 @ X2.T
    return signal_var * np.exp(-0.5 * sqdist / length_scale**2)

# GP posterior via Cholesky decomposition
K_train = rbf_kernel(X_train, X_train) + sigma_n**2 * np.eye(len(X_train))
K_star = rbf_kernel(X_test, X_train)
K_ss = rbf_kernel(X_test, X_test)

L = np.linalg.cholesky(K_train)                       # O(n^3) factorization
alpha_gp = np.linalg.solve(L.T, np.linalg.solve(L, y_train))  # O(n^2) solve
mu_post = K_star @ alpha_gp                            # posterior mean

v = np.linalg.solve(L, K_star.T)
var_post = np.diag(K_ss) - np.sum(v**2, axis=0)       # posterior variance
std_post = np.sqrt(np.maximum(var_post, 0))            # clamp for numerics

GP Posterior Explorer

Click to add observations (max 10)

Kernel:

RBFMatern 3/2Linear

ℓ = 1.00σ² = 0.10Clear0 observations

---

The Indian Buffet Process

The Dirichlet process provides a nonparametric prior for clustering (each observation belongs to exactly one cluster). But what if we want each observation to possess multiple latent features? The Indian Buffet Process (IBP), introduced by Griffiths and Ghahramani (2005), provides a nonparametric prior on binary feature matrices with infinitely many columns.

---

Posterior Consistency and Contraction Rates

The deepest connection between Bayesian nonparametrics and the PAC learning framework lies in the theory of posterior consistency. Just as PAC learning asks “does the learner converge to a good hypothesis as $n \to \infty$?”, posterior consistency asks “does the Bayesian posterior converge to the truth?”

Schwartz’s Theorem

The classical result on posterior consistency is due to Schwartz (1965), who identified the key condition: the prior must assign positive probability to Kullback-Leibler neighborhoods of the truth.

Posterior Contraction Rates

Modern theory (Ghosal, Ghosh & van der Vaart, 2000) goes beyond consistency to rates. The key result establishes that posterior contraction rates are governed by the interplay between prior concentration (how much mass the prior places near the truth) and model complexity (measured by metric entropy).

Consistency of DP Mixture Models

---

Connections & Further Reading

Connection Map

Topic	Domain	Relationship
PAC Learning Framework	Probability & Statistics	Posterior contraction rates parallel PAC sample complexity; Bayesian model selection provides an alternative to SRM
Concentration Inequalities	Probability & Statistics	Posterior contraction proofs use concentration of the log-likelihood ratio; GP concentration bounds use sub-Gaussian tail inequalities
Measure-Theoretic Probability	Probability & Statistics	The DP is defined on measure spaces; posterior consistency uses dominated convergence and the law of large numbers
Spectral Theorem	Linear Algebra	GP kernel matrices are positive semi-definite; the eigendecomposition of the kernel determines the GP’s RKHS
SVD	Linear Algebra	Low-rank GP approximations (Nyström method) use the SVD of the kernel matrix
PCA & Low-Rank Approximation	Linear Algebra	Kernel PCA is equivalent to projecting onto the leading eigenfunctions of the GP kernel; functional PCA uses GP priors

Key Notation Summary

Symbol	Meaning
$\text{DP}(\alpha, G_0)$	Dirichlet process with concentration $\alpha$ and base measure $G_0$
$G_0$	Base measure (prior guess for DP draws)
$\alpha$	Concentration parameter
$\delta_\theta$	Point mass (Dirac delta) at $\theta$
$\text{Dir}(\boldsymbol{\alpha})$	Dirichlet distribution with parameter vector $\boldsymbol{\alpha}$
$V_k \stackrel{\text{iid}}{\sim} \text{Beta}(1, \alpha)$	Stick-breaking beta variables
$w_k = V_k \prod_{j<k}(1-V_j)$	Stick-breaking weight
$\hat{F}_n = \frac{1}{n}\sum_{i=1}^n \delta_{\theta_i}$	Empirical distribution
$K_n$	Number of distinct clusters after $n$ observations
$\mathcal{GP}(m, k)$	Gaussian process with mean $m$ and kernel $k$
$\mathbf{k}_*^\top (\mathbf{K} + \sigma_n^2 \mathbf{I})^{-1} \mathbf{y}$	GP posterior mean
$\text{KL}(p_0 \| p)$	Kullback-Leibler divergence
$N(\varepsilon, \Theta, d)$	$\varepsilon$-covering number of $\Theta$ under metric $d$
$n^{-\beta/(2\beta+d)}$	Minimax contraction rate for $\beta$-smooth densities in $\mathbb{R}^d$