## ----include=FALSE------------------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
library(STARRS)
set.seed(1)

## ----simul, eval=TRUE---------------------------------------------------------
# Dimensions
K <- 3; nk <- rep(500, K); d <- 5
# Gaussian parameters
mu <- rbind(mu1=rep(0, d), mu2=rep(3, d), mu3=rep(-3, d))
sigma <- array(dim=c(K, d, d))
sigma[1, , ] <- 2*diag(d); sigma[2, , ] <- diag(1:d); sigma[3, , ] <- diag(1/(1:d))
# Student contamination
fraction <- 0.3
dfT <- 1; dfT <- 1; muT <- mu; sigmaT <- sigma
# Simulation
GMMcontStudent <- rGMMcontStudent(nk=nk, muG=mu, sigmaG=sigma,
                                  delta=fraction,
                                  dfT=dfT, muT=muT, sigmaT=sigmaT)
print(head(STARRS::GMMcontStudent$X))
print(head(STARRS::GMMcontStudent$Z))
print(head(STARRS::GMMcontStudent$C))

## ----rgmm, eval=TRUE----------------------------------------------------------
robustGMMOutput <- multipleRobustMM(GMMcontStudent$X, nclust=1:5)
print(robustGMMOutput$bestresult$centers)

## ----selection, eval=TRUE-----------------------------------------------------
plot(robustGMMOutput, what='BIC')
plot(robustGMMOutput, what='ICL')

## ----uncertainty, eval=TRUE---------------------------------------------------
plot(robustGMMOutput, what='Two_Dim_Uncertainty')

## ----kmeans, eval=TRUE--------------------------------------------------------
K <- 3; 
# k-means
kmeans <- lapply(1:length(homoGMMStudentCont_K3_d5_n1500), function(f){
  kmeans(homoGMMStudentCont_K3_d5_n1500[[f]]$X, centers=K, nstart=20, iter.max=100)})
# k-medians
kmedians <- lapply(1:length(homoGMMStudentCont_K3_d5_n1500), function(f){
  kmedians(homoGMMStudentCont_K3_d5_n1500[[f]]$X)})

## ----kmeansAri, eval=TRUE-----------------------------------------------------
# Classification accuracy
fractionList <- (0:5)/10; library(mclust)
plot(fractionList, sapply(1:length(fractionList), function(ff){
  adjustedRandIndex(heteroGMMStudentCont_K3_d5_n1500[[ff]]$Z, kmeans[[ff]]$cluster)
  }), type='b', pch=20, xlab='fraction of outliers', ylab='ARI')
points(fractionList, sapply(1:length(fractionList), function(ff){
  adjustedRandIndex(heteroGMMStudentCont_K3_d5_n1500[[ff]]$Z, kmedians[[ff]]$bestresult$cluster)
  }), type='b', pch=20, col='red')

## ----gmm, eval=TRUE-----------------------------------------------------------
# Regular EM
gmm <- lapply(1:length(heteroGMMStudentCont_K3_d5_n1500), function(f){
  Mclust(heteroGMMStudentCont_K3_d5_n1500[[f]]$X, G=K)})
# Robust EM
rgmm <- lapply(1:length(heteroGMMStudentCont_K3_d5_n1500), function(f){
  multipleRobustMM(heteroGMMStudentCont_K3_d5_n1500[[f]]$X, nclust=K)})

## ----gmmAri, eval=TRUE--------------------------------------------------------
fractionList <- (0:5)/10
plot(fractionList, sapply(1:length(fractionList), function(ff){
  adjustedRandIndex(heteroGMMStudentCont_K3_d5_n1500[[ff]]$Z, gmm[[ff]]$classification)
  }), type='b', pch=20, xlab='fraction of outliers', ylab='ARI')
points(fractionList, sapply(1:length(fractionList), function(ff){
  adjustedRandIndex(heteroGMMStudentCont_K3_d5_n1500[[ff]]$Z, rgmm[[ff]]$bestresult$clusters)
  }), type='b', pch=20, col='red')