Dirichlet-Process & Pitman-Yor Mixtures (R)¶
Random measures, bivariate clustering, and the frequentist counterparts¶
The R counterpart to dpmix_python.ipynb. We (1) build a Dirichlet-process draw by stick-breaking from scratch, (2) fit a bivariate DP mixture to Old Faithful with the dirichletprocess package, (3) put it beside its frequentist analogs — mclust (EM finite mixture, chosen by BIC) and a 2-D kernel density estimate (MASS::kde2d) — and (4) contrast the DP with the two-parameter Pitman-Yor process. Data: the faithful geyser (272 eruptions, duration times waiting).
options(repr.plot.width=13, repr.plot.height=4.6)
.libPaths(c("C:/Users/user/R/win-library/4.6", .libPaths()))
suppressMessages({library(dirichletprocess); library(mclust); library(MASS)})
BLUE<-"#2b6cb0"; RED<-"#c53030"; GREEN<-"#2f855a"; ORANGE<-"#dd6b20"; PURP<-"#6b46c1"; GREY<-"#718096"
F <- as.matrix(read.csv("faithful.csv")); colnames(F)<-c("eruptions","waiting")
cat(nrow(F), "eruptions; columns:", paste(colnames(F),collapse=", "), "\n")
Warning message: "package 'dirichletprocess' was built under R version 4.6.1"
272 eruptions; columns: eruptions, waiting
1. Stick-breaking a random measure¶
A draw $G\sim\text{DP}(\alpha,G_0)$ is discrete: $v_k\sim\text{Beta}(1,\alpha)$, $w_k=v_k\prod_{j<k}(1-v_j)$, atoms $\theta_k\sim G_0$. Small $\alpha$ piles the weight on a few atoms; large $\alpha$ spreads it and $G\to G_0$.
stick <- function(alpha, K, d=0){ v<-rbeta(K,1-d,alpha+(1:K)*d); v[K]<-1; v*c(1,cumprod(1-v)[-K]) }
expK <- function(alpha,n,d=0){ ek<-0; for(i in 0:(n-1)) ek<-ek+(alpha+d*ek)/(alpha+i); ek }
set.seed(1)
par(mfrow=c(1,3), mar=c(4,4,3,1))
for(al in c(1,5,25)){ w<-stick(al,40); th<-rnorm(40)
plot(th, w, type="h", col=GREY, lwd=1, xlim=c(-3,3), ylim=c(0,max(w)),
xlab=expression(theta), ylab="weight", main=bquote(alpha==.(al)~" E[K]/200"%~~%.(round(expK(al,200)))))
points(th, w, pch=19, col=RED, cex=.7)
curve(0.06*dnorm(x), add=TRUE, col=GREEN, lwd=1.5) }
par(mfrow=c(1,1))
cat("Small alpha: a few atoms hold nearly all the mass (strong clustering). Large alpha: the weight spreads\n")
cat("across many atoms and the random measure fills in toward the base G0.\n")
Small alpha: a few atoms hold nearly all the mass (strong clustering). Large alpha: the weight spreads
across many atoms and the random measure fills in toward the base G0.
2. The frequentist counterparts -- mclust and a kernel density estimate¶
The classical routes: mclust fits Gaussian finite mixtures by EM and picks the number of components by BIC (the frequentist "how many clusters"); kde2d is the 2-D kernel density estimate (pure smoothing). Both on Old Faithful.
mc <- densityMclust(F, G=1:6, plot=FALSE)
cat("mclust BIC selects G =", mc$G, "components (model", mc$modelName, ")\n")
gx <- seq(1.4,5.4,length=45); gy <- seq(41,98,length=45); grid <- expand.grid(eruptions=gx, waiting=gy)
dmc <- matrix(predict(mc, newdata=grid), 45, 45) # mclust mixture density
kd <- kde2d(F[,1], F[,2], n=45, lims=c(1.4,5.4,41,98)) # kernel density estimate
par(mfrow=c(1,2), mar=c(4,4,3,1))
image(gx,gy,kd$z, col=hcl.colors(20,"viridis"), xlab="eruption (min)", ylab="waiting (min)", main="Kernel density estimate (frequentist)")
points(F, pch=19, cex=.3, col="white")
image(gx,gy,dmc, col=hcl.colors(20,"viridis"), xlab="eruption (min)", ylab="waiting (min)", main=sprintf("mclust EM mixture, G=%d (BIC)", mc$G))
points(F, pch=19, cex=.35, col=c("white","black","red","cyan")[mc$classification])
par(mfrow=c(1,1))
cat("KDE just smooths; mclust also returns a clustering and a component count -- but BIC hands back a single G,\n")
cat("with no distribution over it. That is exactly the gap the Bayesian nonparametric model fills next.\n")
mclust BIC selects G = 3 components (model EEE )
KDE just smooths; mclust also returns a clustering and a component count -- but BIC hands back a single G,
with no distribution over it. That is exactly the gap the Bayesian nonparametric model fills next.
3. The Bayesian nonparametric fit -- dirichletprocess¶
DirichletProcessMvnormal puts a DP prior on a bivariate Normal mixture and samples the clustering, with the number of clusters inferred. We read off its posterior over the number of clusters (the count mclust had to fix) and the clustering itself.
set.seed(2)
dp <- DirichletProcessMvnormal(scale(F)) # standardised inside
dp <- Fit(dp, 800, progressBar=FALSE)
Kchain <- sapply(dp$weightsChain, length) # #clusters at each iteration
Kpost <- Kchain[400:length(Kchain)] # post-burn
cat("dirichletprocess posterior number of clusters: mode", as.integer(names(which.max(table(Kpost)))),
" mean", round(mean(Kpost),2), "\n")
cat(sprintf("(mclust BIC picked %d; this DP puts its mode at %d, and the from-scratch Python sampler puts\n",
mc$G, as.integer(names(which.max(table(Kpost))))))
cat("its mode at 3 -- three engines, three answers, on the same 272 points.)\n")
lab <- dp$clusterLabels; ul <- sort(unique(lab)); cols <- c(BLUE,RED,GREEN,ORANGE,PURP,GREY)
par(mfrow=c(1,2), mar=c(4,4,3,1))
barplot(table(factor(Kpost, levels=min(Kpost):max(Kpost)))/length(Kpost), col=BLUE,
xlab="number of clusters", ylab="posterior", main="DP posterior over #clusters")
plot(F, col=cols[match(lab,ul)], pch=19, cex=.7, xlab="eruption (min)", ylab="waiting (min)", main="DP clustering (representative)")
par(mfrow=c(1,1))
cat("The DP recovers the two main clouds and returns the whole posterior over the number of clusters,\n")
cat("rather than mclust's single BIC-selected value -- the honest count comes with its uncertainty.\n")
cat("\nAnd the disagreement with the other engines is the lesson, not a defect. This package places its\n")
cat("own conjugate base measure on the component means and covariances and samples the concentration\n")
cat("parameter; the from-scratch sampler uses a different normal-inverse-Wishart base and fixes it.\n")
cat("Neither is more correct. The number of clusters a DP mixture reports is a property of the model\n")
cat("and of what you choose to count, not a fact recoverable from the data -- which is precisely why\n")
cat("the POSTERIOR over K is the deliverable and a single selected K is not.\n")
dirichletprocess posterior number of clusters: mode 2 mean 2.36
(mclust BIC picked 3; this DP puts its mode at 2, and the from-scratch Python sampler puts
its mode at 3 -- three engines, three answers, on the same 272 points.)
The DP recovers the two main clouds and returns the whole posterior over the number of clusters,
rather than mclust's single BIC-selected value -- the honest count comes with its uncertainty.
And the disagreement with the other engines is the lesson, not a defect. This package places its
own conjugate base measure on the component means and covariances and samples the concentration
parameter; the from-scratch sampler uses a different normal-inverse-Wishart base and fixes it.
Neither is more correct. The number of clusters a DP mixture reports is a property of the model
and of what you choose to count, not a fact recoverable from the data -- which is precisely why
the POSTERIOR over K is the deliverable and a single selected K is not.
4. Dirichlet process vs Pitman-Yor -- power-law clustering¶
Pitman-Yor adds a discount $d$: join an existing cluster with weight $\propto n_k-d$, open a new one with weight $\propto\alpha+dK$. The number of clusters then grows like a power law $n^d$ instead of the DP's $\log n$ -- the right prior when there are many rare types.
ns <- c(10,25,50,100,200,400,800,1600)
curves <- list("DP d=0"=sapply(ns,function(n) expK(1,n,0)),
"PY d=0.5"=sapply(ns,function(n) expK(1,n,0.5)),
"PY d=0.7"=sapply(ns,function(n) expK(1,n,0.7)))
cols3 <- c(BLUE,RED,ORANGE)
par(mar=c(4,4,3,1))
plot(ns, curves[[1]], type="b", pch=19, log="x", col=cols3[1], ylim=c(0,max(unlist(curves))),
xlab="n (log scale)", ylab="E[ #clusters ]", main="DP grows like log n; Pitman-Yor like a power law")
for(i in 2:3) lines(ns, curves[[i]], type="b", pch=19, col=cols3[i])
legend("topleft", names(curves), col=cols3, lwd=2, pch=19, bty="n")
cat("The discount d makes new clusters cheaper, so Pitman-Yor keeps opening them: its cluster count climbs as a\n")
cat("power of n, the DP only logarithmically. For heavy-tailed, many-rare-types data (words, species) that is the\n")
cat("difference between a prior that keeps up and one that saturates.\n")
The discount d makes new clusters cheaper, so Pitman-Yor keeps opening them: its cluster count climbs as a
power of n, the DP only logarithmically. For heavy-tailed, many-rare-types data (words, species) that is the
difference between a prior that keeps up and one that saturates.
5. Summary¶
Three views of the same idea in R. Stick-breaking shows a Dirichlet process is a discrete random measure with concentration $\alpha$ as its clustering dial. dirichletprocess fits the bivariate DP mixture to Old Faithful and returns a posterior over the number of clusters -- where the frequentist mclust (EM + BIC) and the kernel density estimate deliver a single component count and a smooth surface with no uncertainty attached. Pitman-Yor replaces the DP's logarithmic cluster growth with a power law.
Together with dpmix_python.ipynb (the from-scratch collapsed sampler, the PyMC stick-breaking cross-check, and the KDE comparison) this is the multivariate, random-measure face of the Dirichlet process -- complementing the univariate "how many components" treatment in the variable-selection arc and the discrete DP-LCA in the latent-class arc. The hierarchical Dirichlet process comes next: several groups sharing one set of atoms.