fix many problems (models appearing twice, irrelevant coefficients in a relevant...
[valse.git] / pkg / R / constructionModelesLassoMLE.R
... / ...
CommitLineData
1#' constructionModelesLassoMLE
2#'
3#' TODO: description
4#'
5#' @param ...
6#'
7#' @return ...
8#'
9#' export
10constructionModelesLassoMLE = function(phiInit, rhoInit, piInit, gamInit, mini, maxi,
11 gamma, X, Y, thresh, tau, S, ncores=3, fast=TRUE, verbose=FALSE)
12{
13 if (ncores > 1)
14 {
15 cl = parallel::makeCluster(ncores, outfile='')
16 parallel::clusterExport( cl, envir=environment(),
17 varlist=c("phiInit","rhoInit","gamInit","mini","maxi","gamma","X","Y","thresh",
18 "tau","S","ncores","verbose") )
19 }
20
21 # Individual model computation
22 computeAtLambda <- function(lambda)
23 {
24 if (ncores > 1)
25 require("valse") #nodes start with an empty environment
26
27 if (verbose)
28 print(paste("Computations for lambda=",lambda))
29
30 n = dim(X)[1]
31 p = dim(phiInit)[1]
32 m = dim(phiInit)[2]
33 k = dim(phiInit)[3]
34 sel.lambda = S[[lambda]]$selected
35# col.sel = which(colSums(sel.lambda)!=0) #if boolean matrix
36 col.sel <- which( sapply(sel.lambda,length) > 0 ) #if list of selected vars
37 if (length(col.sel) == 0)
38 return (NULL)
39
40 # lambda == 0 because we compute the EMV: no penalization here
41 res = EMGLLF(phiInit[col.sel,,],rhoInit,piInit,gamInit,mini,maxi,gamma,0,
42 X[,col.sel], Y, tau, fast)
43
44 # Eval dimension from the result + selected
45 phiLambda2 = res$phi
46 rhoLambda = res$rho
47 piLambda = res$pi
48 phiLambda = array(0, dim = c(p,m,k))
49 for (j in seq_along(col.sel))
50 phiLambda[col.sel[j],sel.lambda[[j]],] = phiLambda2[j,sel.lambda[[j]],]
51 dimension = length(unlist(sel.lambda))
52
53 # Computation of the loglikelihood
54 densite = vector("double",n)
55 for (r in 1:k)
56 {
57 if (length(col.sel)==1){
58 delta = (Y%*%rhoLambda[,,r] - (X[, col.sel]%*%t(phiLambda[col.sel,,r])))
59 } else delta = (Y%*%rhoLambda[,,r] - (X[, col.sel]%*%phiLambda[col.sel,,r]))
60 densite = densite + piLambda[r] *
61 det(rhoLambda[,,r])/(sqrt(2*base::pi))^m * exp(-diag(tcrossprod(delta))/2.0)
62 }
63 llhLambda = c( sum(log(densite)), (dimension+m+1)*k-1 )
64 list("phi"= phiLambda, "rho"= rhoLambda, "pi"= piLambda, "llh" = llhLambda)
65 }
66
67 # For each lambda, computation of the parameters
68 out =
69 if (ncores > 1)
70 parLapply(cl, 1:length(S), computeAtLambda)
71 else
72 lapply(1:length(S), computeAtLambda)
73
74 if (ncores > 1)
75 parallel::stopCluster(cl)
76
77 out
78}