projects / valse.git / blame
| Commit | Line | Data |
|---|---|---|
| 2279a641 BA |
1 | #' constructionModelesLassoMLE |
| 2 | #' | |
| 3 | #' TODO: description | |
| 4 | #' | |
| 5 | #' @param ... | |
| 6 | #' | |
| 7 | #' @return ... | |
| 8 | #' | |
| 9 | #' export | |
| 10 | constructionModelesLassoMLE = function(phiInit, rhoInit, piInit, gamInit, mini, maxi, | |
| 11 | gamma, X, Y, seuil, tau, selected, ncores=3, verbose=FALSE) | |
| 46a2e676 | 12 | { |
| 2279a641 BA |
13 | if (ncores > 1) |
| 14 | { | |
| 15 | cl = parallel::makeCluster(ncores) | |
| 16 | parallel::clusterExport( cl, envir=environment(), | |
| 17 | varlist=c("phiInit","rhoInit","gamInit","mini","maxi","gamma","X","Y","seuil", | |
| 18 | "tau","selected","ncores","verbose") ) | |
| 19 | } | |
| 20 | ||
| 21 | # Individual model computation | |
| 22 | computeAtLambda <- function(lambda) | |
| 23 | { | |
| 24 | if (ncores > 1) | |
| 25 | require("valse") #// nodes start with an ampty 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 | ||
| 35 | sel.lambda = selected[[lambda]] | |
| 36 | # col.sel = which(colSums(sel.lambda)!=0) #if boolean matrix | |
| 37 | col.sel <- which( sapply(sel.lambda,length) > 0 ) #if list of selected vars | |
| 38 | ||
| 39 | if (length(col.sel) == 0) | |
| 40 | return (NULL) | |
| 41 | ||
| 42 | # lambda == 0 because we compute the EMV: no penalization here | |
| 43 | res = EMGLLF(phiInit[col.sel,,],rhoInit,piInit,gamInit,mini,maxi,gamma,0, | |
| 44 | X[,col.sel],Y,tau) | |
| 45 | ||
| 46 | # Eval dimension from the result + selected | |
| 47 | phiLambda2 = res_EM$phi | |
| 48 | rhoLambda = res_EM$rho | |
| 49 | piLambda = res_EM$pi | |
| 51485a7d | 50 | phiLambda = array(0, dim = c(p,m,k)) |
| 2279a641 BA |
51 | for (j in seq_along(col.sel)) |
| 52 | phiLambda[col.sel[j],,] = phiLambda2[j,,] | |
| 53 | ||
| 54 | dimension = 0 | |
| 55 | for (j in 1:p) | |
| 56 | { | |
| 57 | b = setdiff(1:m, sel.lambda[,j]) | |
| 58 | if (length(b) > 0) | |
| 59 | phiLambda[j,b,] = 0.0 | |
| 60 | dimension = dimension + sum(sel.lambda[,j]!=0) | |
| 61 | } | |
| 62 | ||
| 63 | # on veut calculer la vraisemblance avec toutes nos estimations | |
| 64 | densite = vector("double",n) | |
| 65 | for (r in 1:k) | |
| 66 | { | |
| 67 | delta = Y%*%rhoLambda[,,r] - (X[, col.sel]%*%phiLambda[col.sel,,r]) | |
| 68 | densite = densite + piLambda[r] * | |
| 69 | det(rhoLambda[,,r])/(sqrt(2*base::pi))^m * exp(-tcrossprod(delta)/2.0) | |
| 70 | } | |
| 71 | llhLambda = c( sum(log(densite)), (dimension+m+1)*k-1 ) | |
| 72 | list("phi"= phiLambda, "rho"= rhoLambda, "pi"= piLambda, "llh" = llhLambda) | |
| 73 | } | |
| 74 | ||
| 75 | #Pour chaque lambda de la grille, on calcule les coefficients | |
| 76 | out = | |
| 77 | if (ncores > 1) | |
| 0eb161e3 | 78 | parLapply(cl, glambda, computeAtLambda) |
| 2279a641 | 79 | else |
| 0eb161e3 | 80 | lapply(glambda, computeAtLambda) |
| 2279a641 BA |
81 | |
| 82 | if (ncores > 1) | |
| 83 | parallel::stopCluster(cl) | |
| 84 | ||
| 85 | out | |
| c3bc4705 | 86 | } |