In Mufabo/ICASSP20.T6.R: ICASSP20.T6.R

knitr::opts_chunk$set(
collapse = TRUE,
comment = "#>"
)

library(ICASSP20.T6.R)
library(ClusterR)

User Input

MC <- 5 # number of Monte Carlo iterations
epsilon <- 0.04 # percantage of replacement outliers
N_k <- 100 # Number of samples per cluster

em_bic <- matrix(c(1,1, 2,2, 2,4, 3,3, 3,4),5, 2, byrow = TRUE)
embic_iter = nrow(em_bic)
nu <- 3 # t
qH <- 0.8 # Huber
cT <- 4.685 # Tukey

data_31

tmp <- data_31(N_k, epsilon)
data <- tmp$data
labels_true <- tmp$labels
r <- tmp$r
N <- tmp$N
K_true <- tmp$K_true
mu_true <- tmp$mu_true
S_true <- tmp$S_true
L_max <- 2 * K_true # search range

Huber Parameters

cH <- sqrt(stats::qchisq(qH, r))
bH <- stats::pchisq(cH^2, r+2) + cH^2/r*(1-stats::pchisq(cH^2, r))
aH <- gamma(r/2)/pi^(r/2) / ( (2*bH)^(r/2)*(gamma(r/2) - pracma::incgam(r/2, cH^2/(2*bH))) + (2*bH*cH^2*exp(-cH^2/(2*bH)))/(cH^2 - bH * r))

Functional Parameters

g <- list(gaus = function(t) g_gaus(t, r),
        t = function(t) g_t(t, r, nu),
      huber = function(t) g_huber(t, r, list(cH, bH, aH)))

rho <- list(gaus = function(t) rho_gaus(t, r),
        t = function(t) rho_t(t, r, nu),
        huber = function(t) rho_huber(t, r, list(cH, bH, aH)),
        tukey = function(t) rho_tukey(t, r, cT)
        )

psi <- list(gaus = function(t) psi_gaus(t),
        t = function(t) psi_t(t, r, nu),
        huber = function(t) psi_huber(t, r, list(cH, bH)),
        tukey = function(t) psi_tukey(t, cT)
        )

eta <- list(gaus = function(t) eta_gaus(t),
        t = function(t) eta_t(t, r, nu),
        huber = function(t) eta_huber(t, r, list(cH, bH)),
        tukey = function(t) eta_tukey(t, cT)
        )

Cluster enumeration

bic <- array(0, c(MC, L_max, 3, embic_iter))
like <- array(0, c(MC, L_max, 3, embic_iter))
pen <- array(0, c(MC, L_max, 3, embic_iter))

oldw <- getOption("warn")
options(warn = -1)

for(iMC in 1:MC){
for(ii_embic in 1:embic_iter){
for(ll in 1:L_max){
 sprintf('iMc: %i, ii_embic: %i, ll: %i', iMC, ii_embic, ll)

 # EM
 tmp <- EM_RES(data, ll, g[[em_bic[ii_embic, 1]]], psi[[em_bic[ii_embic, 1]]])

 mu_est <- tmp$mu_hat
 S_est <- tmp$S_hat
 t <- tmp$t
 R <- tmp$R
 mem <- (R == apply(R, 1, max))

 # BIC
 bicf <- BIC_F(data, S_est, mu_est, t, mem, rho[[em_bic[ii_embic, 2]]], psi[[em_bic[ii_embic, 2]]], eta[[em_bic[ii_embic, 2]]])
 bica <- BIC_A(S_est, t, mem, rho[[em_bic[ii_embic, 2]]], psi[[em_bic[ii_embic, 2]]], eta[[em_bic[ii_embic, 2]]])
 bics <- BIC_S(S_est, t, mem, rho[[em_bic[ii_embic, 2]]])

 bic[iMC, ll, , ii_embic] <- c(bicf$bic, bica$bic, bics$bic)
 like[iMC, ll, , ii_embic] <- c(bicf$like, bica$like, bics$like)
 pen[iMC, ll, , ii_embic] <- c(bicf$pen, bica$pen, bics$pen)
}
} 
}

options(warn = oldw)

Averaging over MC

bic_avg <- rowMeans(aperm(bic, c(2,3,4,1)), dims=3)
pen_avg <- rowMeans(aperm(pen, c(2,3,4,1)), dims=3)
like_avg <- rowMeans(aperm(like, c(2,3,4,1)), dims=3)

Plots

ICASSP20.T6.R::plot_scatter(cbind(labels_true, data), K_true, r)

marker = c('o','s','d','*','x','^','v','>','<','p','h', '+','o')
names = c("Finite", "Asymptotic", "Schwarz")
g_names = c("Gaus", "t", "Huber", "Tukey")

BIC

for(ii_embic in 1:embic_iter){
graphics::matplot(bic_avg[,,ii_embic], lwd = 1.5, xlab = "number of clusters", ylab = "BIC", pch=c("F", "A", "S"), type = 'b', col=1:3)
graphics::title(paste("Nk:",toString(N_k),", eps:", toString(epsilon),", EM-", g_names[em_bic[ii_embic,1]], ", BIC-", g_names[em_bic[ii_embic,2]]))


graphics::grid()
graphics::legend("topleft", legend=names, lty=1:3, cex=0.8, col=1:3)

}

Likelihood

graphics::matplot(like_avg[,1,], lwd = 1.5, xlab = "number of clusters", ylab = "Likelihood", pch=marker[1:5], type = 'b', col=1:5)
graphics::title(paste("Nk:",toString(N_k),", eps:", toString(epsilon)))


leg_names = c()
for (i in 1:embic_iter) {
  leg_names <- c(leg_names, paste("EM-", g_names[em_bic[i, 1]], ", BIC-", g_names[em_bic[i,2]]))
}
graphics::grid()
graphics::legend("topleft", legend=leg_names, lty=1:5, pch=marker[1:5], cex=0.8, col=1:5)

Penalty

for(ii_embic in 1:embic_iter){
graphics::matplot(pen_avg[,,ii_embic], lwd = 1.5, xlab = "number of clusters", ylab = "Penalty", pch=c("F", "A", "S"), type = 'b', col=1:3)
graphics::title(paste("Nk:",toString(N_k),", eps:", toString(epsilon),", EM-", g_names[em_bic[ii_embic,1]], ", BIC-", g_names[em_bic[ii_embic,2]]))


graphics::grid()
graphics::legend("topleft", legend=names, lty=1:3, cex=0.8, col=1:3)

}