## Replication materials for JSS: MixGHD ##

## Figure 1a
library("MixGHD")
load("dataset1.Rdata")
set.seed(12345)
output <- MSGHD(X, G = 3, max.iter = 10000, eps = 0.1, method = "random",
                scale = FALSE)
contourpl(output)

## Figure 1b
load("dataset1.Rdata")
set.seed(12345)
output2 <- cMSGHD(X, G = 3, max.iter = 10000, eps = 0.1, method = "random",
                  scale = FALSE, modelSel = "BIC")
contourpl(output2)

## 7. MixGHD R package

library("MixGHD")
data(bankruptcy, package = "MixGHD")
res <- MCGHD(data = bankruptcy[, 2:3], G = 2:3, max.iter = 1000, method = "kmedoids",
             modelSel = "BIC")
summary(res)

ls(res@gpar[[1]])
ls(res@gpar[[2]])
ARI(res@map, bankruptcy[, 1])
table(res@map, bankruptcy[, 1])


## Figure 2a
plot(bankruptcy[, 2:3], col = res@map, pch = res@map)

## Figure 2b
plot(res@loglik, ylab = "log-likelihood")


## Figure 3a
res1 <- MGHD(data = bankruptcy[, 2:3], G = 2:3, max.iter = 1000, method = "kmedoids",
             modelSel = "BIC")
contourpl(res1)

ARI(res1@map, bankruptcy[, 1])


## Figure 3b
plot(res1)
plot(res)


## 7.1. Data generation and density estimation
set.seed(12345)
data1 <- rCGHD(n = 600, p = 2)

set.seed(12345)
data2 <- rCGHD(n = 600, p = 2, alpha = c(2, -2), omegav = c(2, 2), omega = 3,
               lambdav = c(0.7, 0.9))

densities <- dCGHD(data2[, 1:2], p = 2, alpha = c(2, -2), omegav = c(2, 2), omega = 3, 
                   lambdav = c(0.7, 0.9))
head(densities, n = 3)

## Figure 4
plot(data1, xlab = "Variable 1", ylab = "Variable 2")

plot(data2, xlab = "Variable 1", ylab = "Variable 2")


## 7.2. Discriminant Analysis
data("sonar", package = "MixGHD")
lab <- as.numeric(factor(sonar[, 61]))

## Divide the dataset into a training and a test set, 
## 30% of each class belong to the test set
test <- sonar[c(1:29, 175:33), 1:60]
testL <- lab[c(1:29, 175:33)]
train <- sonar[c(30:174), 1:60]
trainL <- lab[c(30:174)]

set.seed(7)
modelDA <- DA(train, trainL, test, testL, max.iter = 200)
ls(modelDA)
modelDA$ARItest
modelDA$ARItrain

## 7.3. Classification
# install.packages("pgmm")
data("wine", package = "pgmm")
lab <- as.numeric(factor(wine[, 1]))
lab[seq(1, 178, 4)] <- 0

resMGHD <- MGHD(wine[, 2:28], G = 3, label = lab, method = "kmedoids")
resMSGHD <- MSGHD(wine[, 2:28], G = 3, label = lab, method = "kmedoids")
rescMSGHD <- cMSGHD(wine[, 2:28], G = 3, label = lab, method = "kmedoids")
resMCGHD <- MCGHD(wine[, 2:28], G = 3, label = lab, method = "kmedoids")

## ARI for MGHD
ARI(resMGHD@map[lab == 0], wine[lab == 0, 1])
## ARI for MSGHD
ARI(resMSGHD@map[lab == 0], wine[lab == 0, 1])
## ARI for cMSGHD
ARI(rescMSGHD@map[lab == 0], wine[lab == 0, 1])
## ARI for MCGHD
ARI(resMCGHD@map[lab == 0], wine[lab == 0, 1])