Commit | Line | Data |
---|---|---|
228ee602 | 1 | #' constructionModelesLassoMLE |
2 | #' | |
3 | #' Construct a collection of models with the Lasso-MLE procedure. | |
4 | #' | |
5 | #' @param phiInit an initialization for phi, get by initSmallEM.R | |
6 | #' @param rhoInit an initialization for rho, get by initSmallEM.R | |
7 | #' @param piInit an initialization for pi, get by initSmallEM.R | |
8 | #' @param gamInit an initialization for gam, get by initSmallEM.R | |
9 | #' @param mini integer, minimum number of iterations in the EM algorithm, by default = 10 | |
10 | #' @param maxi integer, maximum number of iterations in the EM algorithm, by default = 100 | |
11 | #' @param gamma integer for the power in the penaly, by default = 1 | |
12 | #' @param X matrix of covariates (of size n*p) | |
13 | #' @param Y matrix of responses (of size n*m) | |
14 | #' @param eps real, threshold to say the EM algorithm converges, by default = 1e-4 | |
15 | #' @param S output of selectVariables.R | |
16 | #' @param ncores Number of cores, by default = 3 | |
17 | #' @param fast TRUE to use compiled C code, FALSE for R code only | |
18 | #' @param verbose TRUE to show some execution traces | |
19 | #' | |
20 | #' @return a list with several models, defined by phi, rho, pi, llh | |
21 | #' | |
22 | #' @export | |
23 | constructionModelesLassoMLE <- function(phiInit, rhoInit, piInit, gamInit, mini, | |
24 | maxi, gamma, X, Y, eps, S, ncores = 3, fast, verbose) | |
25 | { | |
26 | if (ncores > 1) | |
27 | { | |
28 | cl <- parallel::makeCluster(ncores, outfile = "") | |
29 | parallel::clusterExport(cl, envir = environment(), varlist = c("phiInit", | |
30 | "rhoInit", "gamInit", "mini", "maxi", "gamma", "X", "Y", "eps", "S", | |
31 | "ncores", "fast", "verbose")) | |
32 | } | |
33 | ||
34 | # Individual model computation | |
35 | computeAtLambda <- function(lambda) | |
36 | { | |
37 | if (ncores > 1) | |
38 | require("valse") #nodes start with an empty environment | |
39 | ||
40 | if (verbose) | |
41 | print(paste("Computations for lambda=", lambda)) | |
42 | ||
9cb34faf BA |
43 | n <- nrow(X) |
44 | p <- ncol(X) | |
45 | m <- ncol(Y) | |
46 | k <- length(piInit) | |
228ee602 | 47 | sel.lambda <- S[[lambda]]$selected |
48 | # col.sel = which(colSums(sel.lambda)!=0) #if boolean matrix | |
49 | col.sel <- which(sapply(sel.lambda, length) > 0) #if list of selected vars | |
50 | if (length(col.sel) == 0) | |
51 | return(NULL) | |
52 | ||
53 | # lambda == 0 because we compute the EMV: no penalization here | |
9cb34faf BA |
54 | res <- EMGLLF(array(phiInit,dim=c(p,m,k))[col.sel, , ], rhoInit, piInit, gamInit, |
55 | mini, maxi, gamma, 0, as.matrix(X[, col.sel]), Y, eps, fast) | |
228ee602 | 56 | |
57 | # Eval dimension from the result + selected | |
58 | phiLambda2 <- res$phi | |
59 | rhoLambda <- res$rho | |
60 | piLambda <- res$pi | |
61 | phiLambda <- array(0, dim = c(p, m, k)) | |
62 | for (j in seq_along(col.sel)) | |
63 | phiLambda[col.sel[j], sel.lambda[[j]], ] <- phiLambda2[j, sel.lambda[[j]], ] | |
64 | dimension <- length(unlist(sel.lambda)) | |
65 | ||
66 | ## Computation of the loglikelihood | |
67 | # Precompute det(rhoLambda[,,r]) for r in 1...k | |
9cb34faf | 68 | detRho <- sapply(1:k, function(r) gdet(rhoLambda[, , r])) |
228ee602 | 69 | sumLogLLH <- 0 |
70 | for (i in 1:n) | |
71 | { | |
72 | # Update gam[,]; use log to avoid numerical problems | |
73 | logGam <- sapply(1:k, function(r) { | |
74 | log(piLambda[r]) + log(detRho[r]) - 0.5 * | |
75 | sum((Y[i, ] %*% rhoLambda[, , r] - X[i, ] %*% phiLambda[, , r])^2) | |
76 | }) | |
77 | ||
78 | logGam <- logGam - max(logGam) #adjust without changing proportions | |
79 | gam <- exp(logGam) | |
228ee602 | 80 | norm_fact <- sum(gam) |
81 | sumLogLLH <- sumLogLLH + log(norm_fact) - log((2 * base::pi)^(m/2)) | |
82 | } | |
83 | llhLambda <- c(sumLogLLH/n, (dimension + m + 1) * k - 1) | |
228ee602 | 84 | list(phi = phiLambda, rho = rhoLambda, pi = piLambda, llh = llhLambda) |
85 | } | |
86 | ||
87 | # For each lambda, computation of the parameters | |
88 | out <- | |
89 | if (ncores > 1) { | |
90 | parLapply(cl, 1:length(S), computeAtLambda) | |
91 | } else { | |
92 | lapply(1:length(S), computeAtLambda) | |
93 | } | |
94 | ||
95 | if (ncores > 1) | |
96 | parallel::stopCluster(cl) | |
97 | ||
98 | out | |
99 | } |