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

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

library(ICASSP20.T6.R)
library(zeallot)
library(parallel)
library(doParallel)
library(foreach)

# User Data
MC <- 5 # number of Monte Carlo iterations
epsilon <- 0.04 # percantage of replacement outliers
N_k <- 50 # 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

out_range <- matrix(c(-20, -20, 20, 20), 2, 2) # range of outliers
step_eps <- 20 #10 # steps between outliers

# Data Generation
x <- seq(out_range[1, 1], out_range[1, 2], step_eps)
y <- seq(out_range[2, 1], out_range[2, 2], step_eps)
c(X, Y) %<-% pracma::meshgrid(x, y)
eps_iter <- length(x)^2

data <- array(0, c(N_k*3,2,eps_iter,MC))
for(ii_eps in 1:eps_iter){
  for(ii_mc in 1:MC){

    c(data[,,ii_eps, ii_mc], labels, r, N, K_true, mu_true, S_true) %<-% data_31(N_k, 0)
    N_repl <- 1
    index_repl <- pracma::randperm(N, N_repl)
    data[index_repl, , ii_eps, ii_mc] <- c(X[ii_eps], Y[ii_eps])
  }
}

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))

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)
        )


bic_final <- array(0, c(MC, eps_iter, L_max, 3, embic_iter))
like_final <- array(0, c(MC, eps_iter, L_max, 3, embic_iter))
pen_final <- array(0, c(MC, eps_iter, L_max, 3, embic_iter))

num_cl <- parallel::detectCores()
cl <- parallel::makeCluster(num_cl)
start_time <- Sys.time()

for(ii_eps in 8:eps_iter){
  # tmp: List(List(bic_1,...),...,List(bic_5,...))
  tmp <- foreach(iMC = 1:MC) %dopar% {
    bic <- array(0, c(L_max, 3, embic_iter))
    like <- array(0, c(L_max, 3, embic_iter))
    pen <- array(0, c(L_max, 3, embic_iter))
    for(iEmBic in 1:embic_iter){
        for(ll in 1:L_max){
          c(mu_est, S_est, t, R) %<-% EM_RES(data[,,ii_eps, iMC], ll, g[[em_bic[iEmBic, 1]]], psi[[em_bic[iEmBic, 1]]])

          mem <- (R == apply(R, 1, max))

          c(bic[ll, 1, iEmBic], pen[ll, 1, iEmBic], like[ll, 1, iEmBic]) %<-% BIC_F(data[,,ii_eps,iMC], S_est, mu_est, t, mem, rho[[em_bic[iEmBic, 2]]], psi[[em_bic[iEmBic, 2]]], eta[[em_bic[iEmBic, 2]]])
          c(bic[ll, 2, iEmBic], pen[ll, 2, iEmBic], like[ll, 2, iEmBic]) %<-% BIC_A(S_est, t, mem, rho[[em_bic[iEmBic, 2]]], psi[[em_bic[iEmBic, 2]]], eta[[em_bic[iEmBic, 2]]])
          c(bic[ll, 3, iEmBic], pen[ll, 3, iEmBic], like[ll, 3, iEmBic]) %<-% BIC_S(S_est, t, mem , rho[[em_bic[iEmBic, 2]]])

        }
    }
    print(c(ii_eps, iMC, Sys.time()-start_time))
    return(list(bic=bic, like=like, pen=pen)) 
  }
  for(i in 1:MC){
    bic_final[i, ii_eps,,,] <- tmp[[i]]$bic
    pen_final[i, ii_eps,,,] <- tmp[[i]]$pen
    like_final[i, ii_eps,,,] <- tmp[[i]]$like
  }

}

Evaluation and Plot

  #####
parallel::stopCluster(cl)
  p_under <- array(0, c(dim(bic_final)[4], eps_iter, embic_iter))
  p_det <- array(0, c(dim(bic_final)[4], eps_iter, embic_iter))
  p_over <- array(0, c(dim(bic_final)[4], eps_iter, embic_iter))

  K_true_det <- pracma::repmat(c(rep(K_true, K_true) == 1:K_true, numeric(L_max-K_true)), MC, 1) == 1

  K_true_under <- pracma::repmat(c(!(rep(K_true, K_true-1) == 1:(K_true-1)), numeric(L_max-K_true+1) ), MC, 1) == 1

  for(iEmBic in 1:embic_iter){
    for(iEpsilon in 1:eps_iter){
      for(k in 1:dim(bic_final)[4]){
        BICmax <- aperm(bic_final[,iEpsilon,,k,iEmBic, drop = FALSE]
                        , c(1, 3, 4, 2, 5)) == apply(aperm(bic_final[,iEpsilon,,k ,iEmBic ,drop = FALSE], c(1,3,4,2,5)), 1, max )

        p_under[k, iEpsilon, iEmBic] <- sum(BICmax[K_true_under])/MC
        p_det[k, iEpsilon, iEmBic] <- sum(BICmax[K_true_det])/MC
        p_over[k, iEpsilon, iEmBic] <- 1 - p_det[k, iEpsilon, iEmBic] - p_under[k, iEpsilon, iEmBic]
      }
    }
  }
#####


names = c("Finite", "Asymptotic", "Schwarz")
g_names = c("Gaus", "t", "Huber", "Tukey")
p_det_2 <- aperm(p_det, c(2,1,3))
c(data2, labels, r, N, K_true, mu_true, S_true) %<-% data_31(N_k, 0)
for(iEmBic in 1: embic_iter){
  for(k_bic in 1:dim(bic_final)[4]){
    layout(t(1:2),widths=c(6,1))
    par(mar=c(4,4,1,0.5))

    plot_scatter(cbind(labels, data2), K_true, r)
    Z <- pracma::Reshape(p_det_2[, k_bic, iEmBic], dim(X)[1], dim(X)[2])
    graphics::contour(x, y, t(Z), add=TRUE, col = grDevices::rainbow(10))
    graphics::title(paste("EM-", g_names[[em_bic[iEmBic, 1]]], ", BIC-", g_names[[em_bic[iEmBic, 2]]], names[k_bic]))

    par(mar=c(.5,1,5,2.5))
    image(z=t(seq(0,1,.1)), col=grDevices::rainbow(10), axes=FALSE, main="Slope", cex.main=.8)
    axis(4,cex.axis=0.8,mgp=c(0,.5,0))
  }
}