################### ## students data ## ################### set.seed(1234) library("flexCWM") data("students", package = "flexCWM") str(students) attach(students) pairs(cbind(WEIGHT, HEIGHT, HEIGHT.F), col = GENDER) fit1 <- cwm(WEIGHT ~ HEIGHT + HEIGHT.F, familyY = list(gaussian, student.t), k = 1:3, initialization = "kmeans") fit2 <- cwm(WEIGHT ~ HEIGHT + HEIGHT.F, familyY = list(gaussian, student.t), Xnorm = cbind(HEIGHT, HEIGHT.F), k = 1:3, initialization = "kmeans", parallel = FALSE) summary(fit2, criterion = "BIC", concomitant = TRUE) plot(fit2) fit2clust <- getCluster(fit2) table(fit2clust, GENDER) ################ ## ExCWM data ## ################ set.seed(17) data("ExCWM", package = "flexCWM") attach(ExCWM) str(ExCWM) resXnorm <- cwm(Xnorm = cbind(Xnorm1, Xnorm2), k = 1:3, initialization = "kmeans") getParXnorm(resXnorm) resXbin <- cwm(Xbin = Xbin, k = 1:3, initialization = "kmeans") getParXbin(resXbin) resXpois <- cwm(Xpois = Xpois, k = 1:3, initialization = "kmeans") getParXpois(resXpois) resX <- cwm(Xnorm = cbind(Xnorm1, Xnorm2), Xbin = Xbin, Xpois = Xpois, k = 1:3, initialization = "kmeans") getParConcomitant(resX) res <- cwm(cbind(Y, Ybtrials-Y) ~ Xnorm1 + Xnorm2 + Xbin + Xpois, data = ExCWM, familyY = binomial, k = 1:3, initialization = "kmeans", Xnorm = cbind(Xnorm1, Xnorm2), Xbin = Xbin, Xpois = Xpois, Xbtrials = Xbtrials[1]) summary(res, concomitant = TRUE) plot(res) ######################### ## Computational times ## ######################### set.seed(1234) library("mnormt") cores <- 1:4 seed <- 10 R <- 100 results <- vector("list", 4) i <- 0 for (n in c(50000, 100000)) { for (d in c(50, 100)) { i <- i + 1 k <- 2 set.seed(seed) results[[i]] <- array(0, c(R, length(cores))) n1 <- n2 <- n/k mu1 <- rep(3, d) mu2 <- -mu1 Sigma1 <- Sigma2 <- diag(d) alpha <- 0 beta1 <- rep(-1, d) beta2 <- rep(1, d) for (r in 1:R) { X1 <- rmnorm(n = n1, mean = mu1, varcov = Sigma1) X2 <- rmnorm(n = n2, mean = mu2, varcov = Sigma2) X <- rbind(X1, X2) Y1 <- rnorm(n/k, mean = (alpha + beta1 %*% t(X[1:(n/k), ])), sd = 0.5) Y2 <- rnorm(n/k, mean = (alpha + beta2 %*% t(X[(n/k+1):n, ])), sd = 0.5) Y <- c(Y1, Y2) loo <- data.frame(Y = Y, X) for (m in cores) { cat(paste0("---> r = ", r, " m=", m, "\n")) options(cl.cores = m) t <- system.time( foo <- cwm(Y ~ ., data = loo, Xnorm = loo[, -1], k = k, initialization = "kmeans", parallel = TRUE, seed = seed) ) results[[i]][r, m] <- t[3] cat(paste0("---> time: ", t[3]), "\n") } } } } b <- ceiling(max(unlist(results))/50)*50 lapply(results, function(res) { boxplot(res, horizontal = TRUE, xlab = "Elapsed time (in seconds)", ylab = "# of used cores", ylim = c(0, b), axes = FALSE) axis(1, at = seq(0, b, by = 50)) axis(2, at = 1:4) })