Number of mixture components — R cross-check¶

Bayesian variable selection, Part 6 (R)¶

Validates dpmix_python.ipynb. The standard R tool for choosing the number of Gaussian mixture components is mclust (Fraley-Raftery): fit mixtures with $G=1,2,\dots$ components by EM and pick $G$ (and the variance structure) by BIC. This is a completely different engine from the Bayesian DP samplers — frequentist EM, not MCMC — so agreement on the number of components is a strong cross-check. It should select about four components, matching the posterior mode.

Data: galaxy_data.csv — 82 galaxy velocities (1000 km/s), Roeder (1990).

In [1]:
.libPaths('C:/Users/user/R/win-library/4.6'); suppressMessages(library(mclust))
y <- read.csv('galaxy_data.csv')$velocity

m <- Mclust(y, G=1:9, verbose=FALSE)                               # EM over 1..9 components, BIC selects G
cat(sprintf('mclust selects G = %d components (variance model "%s")\n', m$G, m$modelName))
bicG <- apply(m$BIC, 1, max, na.rm=TRUE)                           # best BIC at each G
cat('BIC by number of components (higher = better):\n'); print(round(bicG, 1))
cat('\ncross-check -> Python DP mixture: posterior mode K = 4, E[K] ~ 5. Agree.\n')

options(repr.plot.width=11, repr.plot.height=4.2); par(mfrow=c(1,2))
plot(1:9, bicG, type='b', pch=19, col='#2c6fbb', xlab='number of components G', ylab='BIC',
     main='mclust: BIC selects ~4 components'); points(m$G, bicG[m$G], col='firebrick', pch=19, cex=1.6)
d <- densityMclust(y, G=m$G, plot=FALSE)
hist(y, breaks=24, freq=FALSE, col='grey85', border='white', xlab='velocity (1000 km/s)',
     main=sprintf('Fitted %d-component density', m$G))
xs <- seq(min(y)-2, max(y)+2, len=400); lines(xs, predict(d, xs), col='#c0392b', lwd=2.4)

mclust selects G = 4 components (variance model "V")
BIC by number of components (higher = better):
     1      2      3      4      5      6      7      8      9 
-489.5 -462.5 -451.1 -447.0 -459.4 -462.1 -455.7 -464.5 -467.8 

cross-check -> Python DP mixture: posterior mode K = 4, E[K] ~ 5. Agree.
No description has been provided for this image

Same component count, a different engine¶

  • mclust selects four components, with BIC peaking at $G=4$ and 3 and 5–7 close behind — the same mode the Bayesian posterior over $K$ favours (BIC gives a point selection, not a posterior, but the ordering agrees). A frequentist EM+BIC search and two Dirichlet-process MCMC samplers independently agreeing on ~4 components is about as robust as a "how many" answer gets on this famously ambiguous dataset.
  • The fitted density reproduces the multimodal galaxy structure — the dominant supercluster near velocity 20 with smaller groups on either side.

References¶

  • Fraley, C. & Raftery, A. E. (2002). Model-based clustering, discriminant analysis, and density estimation. JASA 97, 611–631.
  • Richardson, S. & Green, P. J. (1997). On Bayesian analysis of mixtures with an unknown number of components. JRSS B 59, 731–792.