How many components? Unknown-dimension mixtures

Bayesian variable selection, Part 6 — Richardson & Green (1997)

The 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.

We use the benchmark galaxy-velocity data (Roeder 1990): the recession velocities of 82 galaxies in the Corona Borealis region. Clusters of galaxies 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 it is exactly the example Congdon (2005) Example 3.8 works through (fitting $K=2,3,4,5$ and comparing). Here we let $K$ be a parameter.

Two standard routes exist for unknown $K$: Green's reversible-jump MCMC (Richardson-Green 1997) with split/merge moves, and the collapsed Dirichlet-process sampler. We build the collapsed one — the same trick as Part 5 — and cross-check against PyMC (a truncated stick-breaking DP) and R (mclust, EM + BIC).

The model and the collapsed 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),\qquad (\mu_k,\sigma_k^2)\sim \text{NIG}(m_0,\kappa_0,a_0,b_0). $$ The Normal-Inverse-Gamma prior is conjugate, so — exactly as the Gamma rates were integrated out in the change-point model — each component's $(\mu_k,\sigma_k^2)$ integrates out analytically, and the predictive density of a point given a component's current members is a closed-form Student-$t$. The sampler then 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}$ (its current size) or opens a brand-new one with probability $\propto \alpha$. As points are reassigned, components are born (a point starts a new one) and die (the last member leaves) — the dimension of the model changes without any Jacobian. This is Neal's (2000) Algorithm 3; it is the collapsed cousin of Green's split/merge RJMCMC.

Hyperparameters (fixed / hyperprior). Base measure $m_0=\bar y$, $\kappa_0=0.05$ (weak on the mean), $a_0=2$, $b_0=\tfrac12\,\mathrm{Var}(y)$. 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: dpmix.py.

import numpy as np, pandas as pd, matplotlib.pyplot as plt
import dpmix                                                        # from-scratch collapsed DP mixture (Neal Alg. 3)

y = pd.read_csv("galaxy_data.csv")["velocity"].values              # 82 galaxy velocities (1000 km/s)
res = dpmix.dp_mixture(y, k0=0.05, a0=2.0, ndraw=10000, burn=5000, seed=1)

kp = res["kpost"]
print("posterior on the number of mixture components K:")
for k in range(3, 8): print("  K=%d : %.3f" % (k, kp[k]))
print("\nE[K] = %.2f   posterior mode = %d   (mean concentration alpha = %.1f)"
      % ((np.arange(len(kp))*kp).sum(), np.argmax(kp), res["alpha"].mean()))

# --- Fig: the data and the posterior-mean mixture density ---
fig, ax = plt.subplots(figsize=(9, 4.4))
ax.hist(y, bins=24, density=True, color='0.8', label='galaxy velocities (n=82)')
ax.plot(res["grid"], res["dens"], color='firebrick', lw=2,
        label='posterior-mean DP-mixture density')
ax.set_xlabel('velocity (1000 km/s)'); ax.set_ylabel('density')
ax.set_title('Galaxy velocities and posterior-mean mixture density')
ax.grid(alpha=.25); ax.legend()
plt.tight_layout(); plt.show()

posterior on the number of mixture components K:
  K=3 : 0.176
  K=4 : 0.278
  K=5 : 0.237
  K=6 : 0.154
  K=7 : 0.085

E[K] = 4.95   posterior mode = 4   (mean concentration alpha = 1.0)

Reading the posterior over dimension

• 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]\approx5$. As with the change-point data, the honest answer is a distribution over the number of components, not a single count; the galaxy data simply do not pin it down to one value.
• The density shows why. The posterior-mean mixture (red) has a dominant central peak near 20 (the main supercluster) with clear shoulders near 9, 16, 23 and a small group near 33 — features that 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 can create or destroy one. The collapsing (integrate out $\mu_k,\sigma_k^2$) is what makes those trans-dimensional moves Jacobian-free.

import pymc as pm, pytensor.tensor as pt                            # cross-check: truncated stick-breaking DP
Kmax = 12; n = len(y)
with pm.Model() as mod:
    alpha = pm.Gamma("alpha", 2, 2)
    beta  = pm.Beta("beta", 1, alpha, shape=Kmax)
    w     = pm.Deterministic("w", beta*pt.concatenate([[1], pt.cumprod(1-beta)[:-1]]))   # stick-breaking weights
    mu    = pm.Normal("mu", y.mean(), 10, shape=Kmax)
    sigma = pm.HalfNormal("sigma", 3, shape=Kmax)
    pm.Mixture("y", w=w, comp_dists=pm.Normal.dist(mu, sigma), observed=y)
    idata = pm.sample(2000, tune=2000, chains=4, target_accept=0.9, random_seed=1, progressbar=False)
W = idata.posterior["w"].values.reshape(-1, Kmax)
K_active = (W > 1.0/n).sum(1)                                       # components holding a non-trivial weight
print("PyMC truncated stick-breaking DP:  E[K_active] = %.1f   mode = %d"
      % (K_active.mean(), np.bincount(K_active).argmax()))
print("-> both samplers agree: about 4 components, with real uncertainty over 3-7.")

# --- Fig: posterior over K, collapsed DP vs stick-breaking DP ---
ks = np.arange(2, 10)
sb = np.bincount(K_active, minlength=int(ks.max()) + 1) / K_active.size
fig, ax = plt.subplots(figsize=(8, 4.4))
width = 0.4
ax.bar(ks - width/2, kp[ks], width, color='firebrick', label='collapsed DP (Neal Alg. 3)')
ax.bar(ks + width/2, sb[ks], width, color='navy', label='truncated stick-breaking DP (PyMC)')
ax.set_xticks(ks)
ax.set_xlabel('number of components K'); ax.set_ylabel('posterior probability')
ax.set_title('Posterior over the number of mixture components')
ax.grid(alpha=.25); ax.legend()
plt.tight_layout(); plt.show()

C:\Users\user\anaconda3\envs\pymc-env\Lib\site-packages\threadpoolctl.py:1226: RuntimeWarning: 
Found Intel OpenMP ('libiomp') and LLVM OpenMP ('libomp') loaded at
the same time. Both libraries are known to be incompatible and this
can cause random crashes or deadlocks on Linux when loaded in the
same Python program.
Using threadpoolctl may cause crashes or deadlocks. For more
information and possible workarounds, please see
    https://github.com/joblib/threadpoolctl/blob/master/multiple_openmp.md

  warnings.warn(msg, RuntimeWarning)
Initializing NUTS using jitter+adapt_diag...

Multiprocess sampling (4 chains in 4 jobs)

NUTS: [alpha, beta, mu, sigma]

Sampling 4 chains for 2_000 tune and 2_000 draw iterations (8_000 + 8_000 draws total) took 24 seconds.

There were 5308 divergences after tuning. Increase `target_accept` or reparameterize.

The rhat statistic is larger than 1.01 for some parameters. This indicates problems during sampling. See https://arxiv.org/abs/1903.08008 for details

The effective sample size per chain is smaller than 100 for some parameters.  A higher number is needed for reliable rhat and ess computation. See https://arxiv.org/abs/1903.08008 for details

PyMC truncated stick-breaking DP:  E[K_active] = 4.5   mode = 4
-> both samplers agree: about 4 components, with real uncertainty over 3-7.

Results — the number of components, inferred not assumed

• Two DP formulations agree. The collapsed Chinese-restaurant sampler and PyMC's truncated stick-breaking both put the mode at $K=4$ with $\mathbb E[K]\approx4.5$–$5$ and mass over 3–7 (the stick-breaking run shows heavy divergences — the multimodal, label-switching geometry mixtures are notorious for — but the component-count summary, which is label-invariant, is stable). And R's mclust (EM + BIC, next notebook) independently selects 4 — three methods, three philosophies, one answer.
• Congdon 3.8 generalised. Where the book fits $K=2,3,4,5$ and compares them, the trans-dimensional sampler returns the whole posterior over $K$ in one run, correctly propagating the uncertainty about how many into the density estimate.
• Closing the unknown-dimension pair. Change-points (Part 5) and mixture components (here) are the two canonical problems where the model's size is unknown. Both are handled by the same recipe — put a conjugate prior on the within-piece parameters, integrate them out, and let a sampler move over the discrete structure (partitions of time, or allocations of points).

References

• Richardson, S. & Green, P. J. (1997). On Bayesian analysis of mixtures with an unknown number of components. JRSS B 59, 731–792.
• Neal, R. M. (2000). Markov chain sampling methods for Dirichlet process mixture models. JCGS 9, 249–265.
• Escobar, M. D. & West, M. (1995). Bayesian density estimation and inference using mixtures. JASA 90, 577–588.
• Roeder, K. (1990). Density estimation with confidence sets exemplified by superclusters and voids in the galaxies. JASA 85, 617–624.
• Congdon, P. (2005). Bayesian Models for Categorical Data, Example 3.8. Wiley.

Data: galaxy_data.csv — velocities of 82 galaxies (1000 km/s), Roeder (1990) / MASS::galaxies. Next: dpmix_R.ipynb — the R cross-check (mclust, EM + BIC over the number of components).