Collection of pivotal algorithms
for: relabelling the MCMC chains in order to undo the label
switching problem in Bayesian mixture models,
as proposed in Egidi, Pappadà, Pauli and Torelli (2018a)
The goal of pivmet
is to propose some pivotal methods in order to:
undo the label switching problem which naturally arises during the MCMC sampling in Bayesian mixture models (\rightarrow) pivotal relabelling (Egidi et al. 2018a)
initialize the K-means algorithm aimed at obtaining a good clustering solution (\rightarrow) pivotal seeding (Egidi et al. 2018b)
You can then install pivmet
from github with:
# install.packages("devtools")devtools::install_github("leoegidi/pivmet")
First of all, we load the package and we import the fish
dataset
belonging to the bayesmix
package:
library(pivmet)#> Loading required package: bayesmix#> Loading required package: runjags#> Loading required package: rjags#> Loading required package: coda#> Linked to JAGS 4.3.0#> Loaded modules: basemod,bugs#> Loading required package: mvtnorm#> Loading required package: RcmdrMisc#> Loading required package: car#> Loading required package: carData#> Loading required package: sandwichdata(fish)y <- fish[,1]N <- length(y) # sample sizek <- 5 # fixed number of clustersnMC <- 12000 # MCMC iterations
Then we fit a Bayesian Gaussian mixture using the piv_MCMC
function:
res <- piv_MCMC(y = y, k = k, nMC = nMC)#> Compiling model graph#> Declaring variables#> Resolving undeclared variables#> Allocating nodes#> Graph information:#> Observed stochastic nodes: 256#> Unobserved stochastic nodes: 268#> Total graph size: 1050#>#> Initializing model
Finally, we can apply pivotal relabelling and inspect the new posterior
estimates with the functions piv_rel
and piv_plot
, respectively:
rel <- piv_rel(mcmc=res, nMC = nMC)piv_plot(y, res, rel, "chains")
piv_plot(y, res, rel, "hist")
Sometimes K-means algorithm does not provide an optimal clustering
solution. Suppose to generate some clustered data and to detect one
pivotal unit for each group with the MUS
(Maxima Units Search
algorithm) function:
#generate some dataset.seed(123)n <- 620centers <- 3n1 <- 20n2 <- 100n3 <- 500x <- matrix(NA, n,2)truegroup <- c( rep(1,n1), rep(2, n2), rep(3, n3))for (i in 1:n1){x[i,]=rmvnorm(1, c(1,5), sigma=diag(2))}for (i in 1:n2){x[n1+i,]=rmvnorm(1, c(4,0), sigma=diag(2))}for (i in 1:n3){x[n1+n2+i,]=rmvnorm(1, c(6,6), sigma=diag(2))}H <- 1000a <- matrix(NA, H, n)for (h in 1:H){a[h,] <- kmeans(x,centers)$cluster}#build the similarity matrixsim_matr <- matrix(1, n,n)for (i in 1:(n-1)){for (j in (i+1):n){sim_matr[i,j] <- sum(a[,i]==a[,j])/Hsim_matr[j,i] <- sim_matr[i,j]}}cl <- KMeans(x, centers)$clustermus_alg <- MUS(C = sim_matr, clusters = cl, prec_par = 5)
Quite often, classical K-means fails in recognizing the true groups:
# launch classical KMeanskmeans_res <- KMeans(x, centers)# plotspar(mfrow=c(1,2))colors_cluster <- c("grey", "darkolivegreen3", "coral")colors_centers <- c("black", "darkgreen", "firebrick")plot(x, col = colors_cluster[truegroup],bg= colors_cluster[truegroup], pch=21,xlab="y[,1]",ylab="y[,2]", cex.lab=1.5,main="True data", cex.main=1.5)plot(x, col = colors_cluster[kmeans_res$cluster],bg=colors_cluster[kmeans_res$cluster], pch=21, xlab="y[,1]",ylab="y[,2]", cex.lab=1.5,main="K-means", cex.main=1.5)points(kmeans_res$centers, col = colors_centers[1:centers],pch = 8, cex = 2)
In such situations, we may need a more robust version of the classical
K-means. The pivots may be used as initial seeds for a classical K-means
algorithm. The function piv_KMeans
works as the classical kmeans
function, with some optional arguments (in the figure below, the colored
triangles represent the pivots).
# launch piv_KMeanspiv_res <- piv_KMeans(x, centers)# plotspar(mfrow=c(1,2), pty="s")colors_cluster <- c("grey", "darkolivegreen3", "coral")colors_centers <- c("black", "darkgreen", "firebrick")plot(x, col = colors_cluster[truegroup],bg= colors_cluster[truegroup], pch=21, xlab="x[,1]",ylab="x[,2]", cex.lab=1.5,main="True data", cex.main=1.5)plot(x, col = colors_cluster[piv_res$cluster],bg=colors_cluster[piv_res$cluster], pch=21, xlab="x[,1]",ylab="x[,2]", cex.lab=1.5,main="piv_Kmeans", cex.main=1.5)points(x[piv_res$pivots[1],1], x[piv_res$pivots[1],2],pch=24, col=colors_centers[1],bg=colors_centers[1],cex=1.5)points(x[piv_res$pivots[2],1], x[piv_res$pivots[2],2],pch=24, col=colors_centers[2], bg=colors_centers[2],cex=1.5)points(x[piv_res$pivots[3],1], x[piv_res$pivots[3],2],pch=24, col=colors_centers[3], bg=colors_centers[3],cex=1.5)points(piv_res$centers, col = colors_centers[1:centers],pch = 8, cex = 2)
Egidi, L., Pappadà, R., Pauli, F. and Torelli, N. (2018a). Relabelling in Bayesian Mixture Models by Pivotal Units. Statistics and Computing, 28(4), 957-969.
Egidi, L., Pappadà, R., Pauli, F., Torelli, N. (2018b). K-means seeding via MUS algorithm. Conference Paper, Book of Short Papers, SIS2018, ISBN: 9788891910233.