Bayesian Nonparametrics
Every model earlier in this collection fixes the size of the question before asking it: how many regression coefficients, how many mixture components, how many knots. Bayesian nonparametrics puts a prior on an infinite-dimensional object instead, and lets the data decide how much of it to use. The name is a misnomer worth clearing up immediately — these models have more parameters than parametric ones, not fewer; what they drop is the commitment to a fixed, finite number of them. The arc runs from priors on clusters (Dirichlet-process and Pitman–Yor mixtures, then hierarchical DPs that share clusters across groups), to priors on functions (Gaussian processes, and their extension to binary outcomes and point processes), to priors on curves (penalised splines and GAMs), on hazards (nonparametric survival), and finally on whole distributions (Pólya trees). Each is built from scratch and set against the frequentist tool that solves the same problem — kernel density estimation, mclust, mgcv, Kaplan–Meier, coxph, ks.test — because the interesting question is rarely which fits better, but what each one is willing to tell you about its own uncertainty.
Two things recur across the seven, and both are worth stating up front. The first is that these models answer a different question than they appear to: the number of clusters a DP mixture reports is a property of the model and the estimand rather than a fact about the data, and three correct engines return three different counts on the same 272 points. The second is that a posterior band prices uncertainty within the assumed structure and is silent about the structure itself — most starkly when the Gaussian process's widening forecast band misses the actual CO₂ level a decade out by nearly six standard deviations, because an RBF kernel extrapolates linearly and the truth is accelerating. Running two independent implementations of every model is what makes both visible — and it is the only reason several of the numbers below can be trusted, since on more than one of these problems two correct-looking engines disagreed and only one of them was right.
What a prior on an infinite-dimensional object is. The phrase sounds forbidding and the idea is not. A prior over coefficients says what values you think a fixed list of numbers might take. A prior over functions says what shapes you think a curve might have — and since a curve is specified by its value at infinitely many points, that is a distribution over an infinite-dimensional space. The same move applies to a distribution (what densities are plausible?), to a hazard (what risk profiles?), and to a clustering (how many groups, and how big?). What makes it workable is that you never handle the infinite object directly: you only ever ask about the finitely many points you actually observed, and each of these priors is constructed so that those finite questions have ordinary, tractable answers. A Gaussian process is the clearest case — any finite set of function values is just a multivariate normal.
Why this is allowed at all. The licence comes from exchangeability: the assumption that the order of your observations carries no information, so that any reshuffling of them is equally probable. De Finetti’s theorem says something remarkable follows — an infinite exchangeable sequence behaves exactly as if the observations were drawn independently from some unknown distribution, with a prior over which distribution that is. The unknown distribution is not an extra assumption you impose; it is forced on you by exchangeability alone. That is why a prior over distributions is the natural object rather than an exotic one, and it is the foundation the Dirichlet process, the Pólya tree and the rest are built on. It also marks the limit: when order does matter — time series, spatial data — exchangeability fails, and the structure has to be put back by hand, which is precisely what a Gaussian process kernel does.
A kernel, and what the marginal likelihood does with it. For the function half of the arc, everything is set by the kernel: a function giving the covariance between the curve’s values at any two inputs. Points close together get high covariance, so the curve varies smoothly; how quickly that covariance decays with distance is the length-scale, and it is the single number controlling wiggliness. Kernels add and multiply, which is what lets the CO₂ model be written as trend plus seasonal plus noise and then read back apart. The reason no cross-validation appears anywhere in these pages is the marginal likelihood — the probability the model assigns to the observed data averaged over every function its prior allows. It penalises a length-scale that is too short (the prior wastes its mass on wiggles the data never show) and one that is too long (the prior cannot reach the data at all), so it selects the hyperparameters on its own. What it cannot do is question the kernel itself, which is why the CO₂ forecast below misses by nearly six standard deviations while its own diagnostics look healthy.
Reading the numbers this arc reports. Four scales recur and none is self-explanatory. The adjusted Rand index compares two clusterings by counting pairs of points the two agree to put together or apart, corrected so that chance agreement scores 0 and a perfect match scores 1 — so an ARI of 0.998 is essentially exact recovery, and the gap between 0.960 and 0.958 is nothing. A log Bayes factor is the log of the ratio of two models’ marginal likelihoods: positive favours the first, and because it is a log, anything past about 5 is already decisive evidence — which is the measure of what a value of 257 means. L₁ error is the area between an estimated density and the truth; it runs from 0 to 2, so 0.19 against 0.31 is a real but not enormous margin. And a divergence is a specific failure of Hamiltonian samplers, where the trajectory destabilises in a region of sharp curvature; divergences cluster in exactly the part of the posterior being missed, which is why 6.2% of draws diverging matters more than the percentage suggests.
The prior that grows, and the prior that cannot
Two panels for the two halves of the arc. A is the stick-breaking construction that Parts 1 and 2 rest on, drawn at the three concentrations the notebook uses. B uses the observed Mauna Loa series together with the 2012 forecasts the two engines publish.
A is the whole idea of an infinite prior in one picture. Break a unit stick at a random point, keep the piece, break what remains, repeat — and you have infinitely many weights that nonetheless sum to one, so a model can carry unlimited components while using only a few. The single parameter α sets how fast they decay, and that decay is the clustering prior: at α = 1 the first three atoms hold 87% of the mass and the data will be pushed into a handful of groups; at α = 25 they hold 11% and the prior is happy to keep opening new ones. Nothing here is truncated or approximated — the weights genuinely continue forever, which is what the word nonparametric is pointing at.
B is the counterweight, and the reason the arc keeps insisting on a second engine.
Both models fit 521 months of CO₂ essentially perfectly and both produce a forecast band that
widens with distance, exactly as it should. Eleven years out only one of them contains the answer.
The Gaussian process says 385.9 ppm, 95% [383.3, 388.6], against an actual
393.8 — outside its own interval by 5.8 standard
deviations. The mgcv GAM on the same series says 391.4 [387.6, 395.2] and
covers, 1.3 standard errors low.
The difference is not tuning and not sample size; it is what each prior does once the data run out. An RBF kernel is stationary, so beyond the last observation it reverts toward the mean and the trend it continues is effectively linear — while atmospheric CO₂ is accelerating. A thin-plate spline extrapolates its polynomial null space instead, which happens to be closer to the truth here. Neither band is wrong about the uncertainty it was asked to price: both correctly report growing ignorance within their own structure, and neither has any way to express doubt about that structure. That distinction is the single most useful thing in this section, and it is why the marginal likelihood choosing the hyperparameters is not the same as the model being right.
How the seven examples relate
Grouped by what the prior is placed on. Each group has a trunk that establishes the construction and a branch that extends it; read the trunk first.
A prior on clusters
How many groups are there, and do separate populations share them?
1 · DP & Pitman–Yor mixturesstick-breaking, and what the concentration and discount buy 2 · Hierarchical DPtwo DPs stacked, so groups draw from one shared menu of clustersA prior on functions
What shapes might the curve have — and what if the outcome is not continuous?
3 · Gaussian-process regressionthe kernel, the marginal likelihood, and the forecast that misses 4 · GP classification & Cox processesthe same latent GP through a link: logit for classes, log for event ratesA prior on everything else
Curves cheaply, hazards without a shape, and whole distributions.
5 · Penalised splines & GAMsa roughness penalty read as a random-walk prior — the cheap route to a GP 6 · Nonparametric survivala free hazard per interval, tied by a Gamma process rather than a formula 7 · Pólya treesa prior over distributions, which doubles as a goodness-of-fit test
Every one is set against the frequentist tool that solves the same problem — a kernel density
estimate, mclust, mgcv, Kaplan–Meier, coxph,
ks.test. The comparison is rarely about which fits better. It is about what each is
willing to say about its own uncertainty, and Part 3 is where that question bites hardest.
Dirichlet-Process & Pitman–Yor Mixtures
The entry point: a prior over discrete random probability measures, built by
stick-breaking and then applied to the Old Faithful eruptions as a bivariate
mixture with the number of components unknown. Set against a kernel density estimate, the
contrast is sharp — the KDE returns a smooth surface and nothing else, while the DP mixture
also delivers a clustering, a posterior over the number of groups, and the
uncertainty in it. Pitman–Yor then adds a discount parameter that makes new
clusters cheaper, turning the DP's logarithmic growth in cluster count into a power law — the
difference between a prior that keeps up with heavy-tailed data and one that saturates.
The project's most useful lesson is one it makes explicit. Three
engines report three different cluster counts on the same 272 points — 3 from
the collapsed Gibbs sampler, 2 from R's dirichletprocess, 4 from PyMC's truncated
stick-breaking — — and the disagreement is not error but definition. They are
not measuring the same thing: one counts occupied clusters, another counts sticks above a weight
threshold, and the PyMC model shares a single covariance across components, which structurally
inflates the count. Each is a correct answer to a slightly different question. The concentration parameter, the obvious suspect, is ruled out by direct test — the
mode holds at 3 across a fourfold range of α. The number of clusters is a property of
the model and the estimand, not a fact about the data, which is exactly why the
posterior over K is the deliverable and a single selected K is not.
(A companion project on galaxy velocities approaches
the same method from the reversible-jump side, in the variable-selection arc.)
The Hierarchical Dirichlet Process
One DP handles one population; real data arrives in groups, and the question is
whether they draw on the same components in different proportions. Separate DPs per group cannot
say — each invents its own labels. Pooling answers by destroying the grouping. The HDP
stacks two DPs: a global one draws a shared menu, each group draws its own
mixing weights over it. Verified against ground truth first, then applied to 342 Palmer penguins
across three islands, where it recovers the species (ARI 0.998) and the
biogeography at once — one cluster on all three islands (Adélie), Gentoo on Biscoe alone,
Chinstrap on Dream alone. The per-group weight table is the whole point of an HDP, and it
reads across groups rather than across clusters — each row is one group's mixing
proportions and sums to one, which is the only normalisation under which the numbers mean
anything. On the simulated data those rows reproduce the generating weights exactly, structural
zeros included. And the honest comparison is worth stating: pooled mclust scores
marginally better on ARI (0.960 against 0.958), so the case for the HDP is not that it
clusters better but that it answers a question the alternatives cannot pose.
Gaussian-Process Regression
From priors on clusters to a prior on functions. Any finite set of function
values is jointly Gaussian with covariance set by a kernel, and with Gaussian
noise the posterior and the marginal likelihood are closed form — no cross-validation appears
anywhere; the marginal likelihood picks the length-scale, amplitude and noise itself. On the
motorcycle data the GP recovers the dip-and-rebound and turns out to be the smoothing spline
in disguise (a spline is a GP posterior mean, which R's mgcv makes
concrete). On Mauna Loa CO₂ the kernel is composed — long-run trend plus a damped
one-year periodic term plus short wiggle — and the decomposition recovers a 56-year trend and
a 6.6 ppm seasonal cycle, matching an independent GAM to a tenth of a ppm.
The finding is in the forecast, and it is a warning. Extrapolated
ten years, the band widens reassuringly — and then fails to cover: the GP says 385.9 ppm
with a 95% interval of [383.3, 388.6] against an actual 393.8, a 5.8σ miss.
A widening band is not the same as an honest one, and only checking it against what happened
tells them apart. The cause is exact rather than vague: an RBF trend reverts to its mean and
extrapolates linearly, and a linear OLS on the same window is 11.5 ppm low while a quadratic
is within 1.4. The band prices uncertainty within the kernel, not about it —
marginal likelihood chooses hyperparameters, not assumptions.
GP Classification & Log-Gaussian Cox Processes
The same latent GP, pushed through a link: logit gives a classifier, log gives
the intensity of a Cox process. The closed form is gone, so the posterior is
handled by a from-scratch Laplace approximation checked against HMC. Validated
first on a known non-monotone probability, where the GP scores RMSE 0.062 against logistic
regression's 0.212 — logistic has log-odds linear in x and is structurally confined to one
monotone S. On the Pima data the notebooks report an honest draw instead of manufacturing a
win: accuracy 0.776 versus 0.773, because the problem is nearly linearly separable. What the
GP adds there is calibrated ignorance — uncertainty that grows in the sparse corners, where
logistic regression stays equally confident having seen almost nobody. The finding is in the Cox process. Fitted in R the coal-disaster intensity came out at 1.5/yr against
Python's 3.0 — a factor-of-two exposure error, and one that looks
perfectly plausible from inside either engine. An offset passed through gam()'s
argument is not carried into predict.gam(), so the returned linear
predictor is already a per-year rate and dividing by the bin width again halves it. The observed
raw rate of 3.15/yr settles which figure was right, and the R notebook now prints the check:
exp(η̂) = 3.055 per year against fitted() = 6.109 per two-year bin, a ratio
of exactly 2.
Bayesian Penalised Splines & Additive Models
The cheap route to the same place a GP reaches: lay down deliberately more B-spline bumps than the data can support, then penalise the roughness of the coefficients. Read as a prior, that penalty is a random walk, which makes the construction Bayesian without changing the algebra — and the smoothing variance is inferred rather than tuned. Against cross-validation the comparison is worth making numerically rather than calling the two "comparable": GCV picks λ = 38, the posterior implies 29, and the two fitted curves differ by 1.3% of the response range. The gap in λ matters so little because the penalty enters only as a ratio inside a matrix inverse. Stacking one spline per covariate gives an additive model whose three partial effects on ozone are physically legible — heat raises it, wind disperses it, sunlight raises then flattens. The sharper observation is about samplers: the PyMC version reproduces the fit to 1.0% of range but 6.2% of its draws diverge, the classic funnel of a hierarchical variance, which the conjugate Gibbs sampler never feels because it draws from exact full conditionals rather than following a gradient through the joint density.
View example →Nonparametric Bayesian Survival
Parametric survival makes you commit to a hazard shape before you look. This one
refuses: a free hazard on every interval, tied together by a Gamma process or
a random walk rather than a formula. Validated first against a known rising hazard
(correlation 0.979), then applied to the Gehan leukemia trial, where the Bayesian survival
curve sits on the Kaplan–Meier step function and the proportional-hazards term recovers the
1963 result — HR 0.19 [0.08, 0.39], 6-MP cutting relapse risk to a fifth of
placebo — while also returning the baseline hazard a Cox fit conditions away. Two cross-engine gaps are worth naming rather than waving through. Python
reports a Cox HR of 0.22 and R 0.21, which is a convention rather than a
disagreement: 18 of the 42 event times are tied, and statsmodels defaults to
the Breslow approximation where coxph defaults to Efron. Fitting both reproduces each
engine's number exactly. The PyMC check is a real gap and reported as one — 0.24 against
0.19, with the baseline hazards 55% apart — which the overlapping intervals accommodate and
30 events make unsurprising.
Pólya Trees
The arc closes on a prior over whole distributions. Split the real line in half, split each half, and put a Beta prior on going left at every node; centre those Betas on a parametric family and one parameter controls how much you trust it. That buys two things at once. As a density estimator it beats a kernel estimate on a bimodal mixture (L₁ error 0.191 against 0.309) and returns a credible band. And because it is conjugate, the marginal likelihood is closed form — so the same object becomes a goodness-of-fit Bayes factor against the family at its centre, agreeing with KS and Shapiro on simulated data but doing what they cannot: handing back the density that fits instead of merely rejecting one. On 3,772 S&P 500 daily returns the log Bayes factor against normality is 257 — the leptokurtosis the Risk and Asset Allocation arc spends eight projects modelling, confirmed here without assuming a Student-t, a GARCH or a copula. The classical tests carry an assumption the Pólya tree does not: KS and Shapiro require a continuous distribution, and 11 of these 3,772 returns are exact zeros. At p = 5e−33 the ties change nothing here, but the Pólya tree never needed the assumption in the first place.
View example →