## ----setup, include=FALSE-----------------------------------------------------
knitr::opts_chunk$set(echo = TRUE)

## ----fig.width=7, fig.height=7------------------------------------------------
library("telescope")
data("diabetes", package = "mclust")
y <- diabetes[, c("glucose", "insulin", "sspg")]
z <- diabetes[, "class"]

pairs(y, col = c("darkred", "blue", "darkgreen")[z], pch = 19)

## -----------------------------------------------------------------------------
priorOnK <- priorOnK_spec("BNB_143")
priorOnWeights <- priorOnAlpha_spec("alpha_const", alpha = 0.5)

## -----------------------------------------------------------------------------
r <- ncol(y)
R <- apply(y, 2, function(x) diff(range(x)))
b0 <- apply(y, 2, median)
B_0 <- rep(1, r)  
B0 <- diag((R^2) * B_0)
c0 <- 2.5 + (r-1)/2
g0 <- 0.5 + (r-1)/2
G0 <- 100 * g0/c0 * diag((1/R^2), nrow = r)
C0 <- g0 * chol2inv(chol(G0))

## -----------------------------------------------------------------------------
M <- 1000
burnin <- 1000
thin <- 1

## -----------------------------------------------------------------------------
set.seed(4711)
z0 <- kmeans(y, centers = 3, nstart = 100)$cluster
mu <- sapply(split(y, z0), colMeans)
Sigma <- array(sapply(split(y, z0), var), c(r, r, ncol(mu)))
eta <- rep(1/ncol(mu), ncol(mu))

## -----------------------------------------------------------------------------
Kmax <- 50
samples <-
    sampleMultNormMixture(
        y, z0, mu, Sigma, eta,
        c0, g0, G0, C0, b0, B0,
        M, burnin, thin, 
        Kmax = Kmax, G = "MixDynamic",
        priorOnK, priorOnWeights = priorOnWeights)

## -----------------------------------------------------------------------------
Mu <- samples$Mu
Eta <- samples$Eta
S <- samples$S
Nk <- samples$Nk
K <- samples$K
Kplus <- samples$Kplus
nonnormpost_mode_list <- samples$nonnormpost_mode

## ----fig.width=7, fig.height=5------------------------------------------------
par(mfrow = c(1, 2))
with(samples,
     matplot(burnin + 1:M, cbind(K, Kplus), type = "l", lty = 1,
             col = c("grey", "black"),
             xlab = "iteration", 
             ylab = expression(`K/` ~K["+"]),
             ylim = c(0, Kmax)))
matplot(burnin + 1:M, samples$Eta, type = "l", lty = 1, col = 1,
        xlab = "iteration", ylim = 0:1,
        ylab = expression(eta["k"]))

## ----fig.width=7, fig.height=7------------------------------------------------
Mu_ <- do.call("rbind", lapply(1:Kmax, function(k) samples$Mu[,,k])) |>
    as.data.frame() |> na.omit()
colnames(Mu_) <- colnames(y)
Mu_ <- subset(Mu_,
              (glucose >= 50 & glucose <= 320) &
              (sspg >= 0 & sspg <= 500) &
              (insulin >= 300 & insulin <= 1300))
pairs(Mu_, col = rgb(0, 0, 0, 0.2), pch = 19)

## -----------------------------------------------------------------------------
(p_Kplus <- tabulate(Kplus, nbins = max(Kplus))/M)

## -----------------------------------------------------------------------------
Kplus_hat <- which.max(p_Kplus)
Kplus_hat

## -----------------------------------------------------------------------------
M0 <- sum(Kplus == Kplus_hat)
M0 

## -----------------------------------------------------------------------------
index <- Kplus == Kplus_hat

Nk[is.na(Nk)] <- 0
Nk_Kplus <- (Nk[index, ] > 0)

Mu_inter <- Mu[index, , , drop = FALSE]
Mu_Kplus <- array(0, dim = c(M0, r, Kplus_hat)) 
for (j in 1:r) {
  Mu_Kplus[, j, ] <- Mu_inter[, j, ][Nk_Kplus]
}

Eta_inter <- Eta[index, ]
Eta_Kplus <- matrix(Eta_inter[Nk_Kplus], ncol = Kplus_hat)

w <- which(index)
S_Kplus <- matrix(0, M0, nrow(y))
for (i in seq_along(w)) {
  m <- w[i]
  perm_S <- rep(0, Kmax)
  perm_S[Nk[m, ] != 0] <- 1:Kplus_hat
  S_Kplus[i, ] <- perm_S[S[m, ]]
}

## -----------------------------------------------------------------------------
Func_init <- t(nonnormpost_mode_list[[Kplus_hat]]$mu)
identified_Kplus <- identifyMixture(
    Func = Mu_Kplus, Mu_Kplus, Eta_Kplus, S_Kplus, Func_init)

## -----------------------------------------------------------------------------
identified_Kplus$non_perm_rate

## ----fig.width=7, fig.height=7------------------------------------------------
Mu_ <- do.call("rbind", lapply(1:Kplus_hat, function(k) identified_Kplus$Mu[,,k])) |>
    as.data.frame()
colnames(Mu_) <- colnames(y)
z_ <- identified_Kplus$class
COLS <- apply(rbind(col2rgb(c("darkred", "blue", "darkgreen")),
                    alpha = 0.2 * 255) / 255,
              2, function(x) do.call("rgb", as.list(x)))
pairs(Mu_, col = COLS[z_], pch = 19)

## ----fig.width=7, fig.height=7------------------------------------------------
colMeans(identified_Kplus$Mu)
colMeans(identified_Kplus$Eta)

z_sp <- apply(identified_Kplus$S, 2,
              function(x) which.max(tabulate(x, Kplus_hat)))
table(z_sp)
table(z, z_sp)

library("mclust")
1 - classError(z_sp, z)$errorRate
adjustedRandIndex(z, z_sp)

pairs(y, col = c("darkred", "blue", "darkgreen")[z_sp], pch = 19)