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

## ----prep, message=FALSE, warning=FALSE,eval = TRUE---------------------------
library(MultiEFM)
library(irlba)   


trace_statistic_fun <- function(H, H0){
  tr_fun <- function(x) sum(diag(x))
  mat1 <- t(H0) %*% H %*% qr.solve(t(H) %*% H) %*% t(H) %*% H0
  return(tr_fun(mat1) / tr_fun(t(H0) %*% H0))
}
trace_list_fun <- function(Hlist, H0list){
  trvec <- sapply(seq_along(Hlist), function(i) trace_statistic_fun(Hlist[[i]], H0list[[i]]))
  return(mean(trvec, na.rm = TRUE))
}

## ----gene, message=FALSE, warning=FALSE,eval = TRUE---------------------------
set.seed(1)
nu <- 2 # nu is set to 2 for heavier tails
p <- 300 # High-dimensional setting
nvec <- c(80, 100);  q <- 3; qs <- c(2,2); S <- length(nvec)
sigma2_eps <- 1
datList <- gendata_simu_robust(seed=1, nvec=nvec, p=p, q=q, qs=qs, rho=c(5,5), err.type='mvt', nu=nu)
XList <- datList$Xlist

## ----fit, message=FALSE, warning=FALSE,eval = TRUE----------------------------
tic <- proc.time()
res <- MultiEFM(XList, q=q, qs_vec=qs, verbose = FALSE)
toc <- proc.time()
time_use <- toc[3] - tic[3]

## ----perf, message=FALSE, warning=FALSE,eval =TRUE----------------------------
results <- data.frame(
  Metric = c('Shared Loading (A)', 'Specific Loading (B)', 'Shared Factor (F)', 'Specific Factor (H)', 'Time (s)'),
  Value = c(trace_statistic_fun(res$A, datList$A0),
            trace_list_fun(res$B, datList$Blist0),
            trace_list_fun(res$F, datList$Flist),
            trace_list_fun(res$H, datList$Hlist),
            time_use)
)
print(results)

## ----auto, message=FALSE, warning=FALSE,eval = TRUE---------------------------
hq_res <- selectFac.MultiEFM(XList, q_max=10, qs_max=5, verbose = FALSE)
message("Estimated shared q = ", hq_res$hq, " VS true q = ", q)

## ----setdown, message=FALSE, warning=FALSE,eval = TRUE------------------------
sessionInfo()