Dirichlet-Process Mixtures — How Many Components?

Python · R · Part 6 of the variable-selection arc — unknown model dimension

Overview

Part 6 of the variable-selection arc: how many components? The reversible-jump change-point sampler of Part 5 answered "how many regimes?"; its sibling asks "how many components generated this data?" — a mixture model whose number of components $K$ is itself unknown. This is the second classic unknown-dimension problem. The benchmark is the galaxy-velocity data (Roeder 1990): the recession velocities of 82 galaxies in Corona Borealis, where clusters at different distances show up as modes in the velocity distribution — so the number of components is a question about the large-scale structure of the universe, and exactly the example Congdon (2005) works through by fitting $K=2,3,4,5$ and comparing. Here we let $K$ be a parameter. We build the collapsed Dirichlet-process sampler — the same integrate-out trick as Part 5 — and cross-check against PyMC (truncated stick-breaking) and R (mclust, EM + BIC).

y_i\mid z_i=k \sim \mathcal N(\mu_k,\sigma_k^2), \quad (\mu_k,\sigma_k^2)\sim\text{NIG}; \qquad \Pr(\text{new component}) \propto \alpha, \quad \Pr(\text{join } k) \propto n_{-i,k}

The collapsed Dirichlet-process sampler

Each observation comes from one of an unknown number of normal components, $y_i\mid z_i=k \sim \mathcal N(\mu_k,\sigma_k^2)$ , with a conjugate Normal-Inverse-Gamma prior on $(\mu_k,\sigma_k^2)$ . Exactly as the Gamma rates were integrated out in the change-point model, each component's parameters integrate out analytically, so the predictive density of a point given a component's current members is a closed-form Student- $t$ and the sampler works only with the discrete allocation $z$ of points to components. A Dirichlet-process (Chinese-restaurant) prior governs the allocation and, crucially, the number of occupied components: a point joins an existing component with probability $\propto n_{-i,k}$ or opens a brand-new one with probability $\propto \alpha$ . As points are reassigned, components are born and die — the dimension of the model changes with no Jacobian. This is Neal's (2000) Algorithm 3, the collapsed cousin of Green's split/merge RJMCMC. The concentration $\alpha$ , which sets the prior on $K$ ( $\mathbb E[K]\approx\alpha\log n$ ), is not fixed but given a $\text{Gamma}(2,2)$ hyperprior and sampled (Escobar–West 1995), so the posterior on $K$ is driven by the data rather than pinned by a chosen $\alpha$ . Engine: a self-contained dpmix.py.

Results

About four components, but genuinely uncertain. The posterior on $K$ peaks at 4 (0.28) with substantial mass on 3, 5 and 6 ( $\mathbb E[K]\approx 5$ ) — as with the change-point data, the honest answer is a distribution over the number of components, not a single count. The posterior-mean density shows why: a dominant central peak near velocity 20 (the main supercluster) with clear shoulders near 9, 16, 23 and a small group near 33 — features three or seven components could also render, hence the spread. This is the mixture face of unknown dimension: where Part 5 moved between numbers of change-points, here the sampler moves between numbers of components, every reassignment able to create or destroy one, and the collapsing is what makes those trans-dimensional moves Jacobian-free. Three methods, three philosophies, one answer: the collapsed Chinese-restaurant sampler and PyMC's truncated stick-breaking both put the mode at $K=4$ ( $\mathbb E[K]\approx 4.5\text{–}5$ , mass over 3–7), and R's mclust (EM + BIC) independently selects 4. Where Congdon fits $K=2,3,4,5$ and compares, the trans-dimensional sampler returns the whole posterior over $K$ in one run, correctly propagating the uncertainty about how many into the density estimate. Together with Part 5 (change-points) and Part 7 (spline knots), this is one of the unknown-dimension trilogy — all handled by the same recipe: conjugate prior on the within-piece parameters, integrate out, and let the sampler move over the discrete structure.

Notebooks

Python Collapsed Dirichlet-process mixture from scratch — Neal Algorithm 3 over the galaxy velocities, inferring the number of components with a sampled concentration, plus a PyMC stick-breaking cross-check (Python) Open → R R cross-check — mclust (EM + BIC) selecting the number of Gaussian components (R) Open →

Downloads

dpmix.py The from-scratch collapsed DP-mixture sampler — NIG-marginalised Student-t predictives, Chinese-restaurant allocation, and Escobar–West concentration sampling (new self-contained module) galaxy_data.csv Recession velocities of 82 galaxies in Corona Borealis (Roeder 1990; MASS::galaxies)

References

Neal, R. M. (2000). Markov chain sampling methods for Dirichlet process mixture models. Journal of Computational and Graphical Statistics 9(2), 249–265. — Algorithm 3, the collapsed sampler reproduced here
Escobar, M. D. & West, M. (1995). Bayesian density estimation and inference using mixtures. Journal of the American Statistical Association 90(430), 577–588. — the concentration-parameter sampling and the galaxy density estimate
Ferguson, T. S. (1973). A Bayesian analysis of some nonparametric problems. Annals of Statistics 1(2), 209–230. — the Dirichlet process
Antoniak, C. E. (1974). Mixtures of Dirichlet processes with applications to Bayesian nonparametric problems. Annals of Statistics 2(6), 1152–1174. — Dirichlet-process mixtures and the prior on the number of components
Richardson, S. & Green, P. J. (1997). On Bayesian analysis of mixtures with an unknown number of components. Journal of the Royal Statistical Society B 59(4), 731–792. — the reversible-jump mixture analysis of the galaxy data
Roeder, K. (1990). Density estimation with confidence sets exemplified by superclusters and voids in the galaxies. Journal of the American Statistical Association 85(411), 617–624. — the galaxy-velocity data
Fraley, C. & Raftery, A. E. (2002). Model-based clustering, discriminant analysis, and density estimation. Journal of the American Statistical Association 97(458), 611–631. — the R mclust EM + BIC cross-check