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, | |
| 08f4604c | 11 | gamma, X, Y, thresh, tau, S, ncores=3, artefact = 1e3, verbose=FALSE) |
| 46a2e676 | 12 | { |
| 08f4604c 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","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 | ||
| 35 | sel.lambda = S[[lambda]]$selected | |
| 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$phi | |
| 48 | rhoLambda = res$rho | |
| 49 | piLambda = res$pi | |
| 50 | phiLambda = array(0, dim = c(p,m,k)) | |
| 51 | for (j in seq_along(col.sel)) | |
| 52 | phiLambda[col.sel[j],,] = phiLambda2[j,,] | |
| 53 | dimension = length(unlist(sel.lambda)) | |
| 54 | ||
| 55 | # Computation of the loglikelihood | |
| 56 | densite = vector("double",n) | |
| 57 | for (r in 1:k) | |
| 58 | { | |
| 59 | delta = (Y%*%rhoLambda[,,r] - (X[, col.sel]%*%phiLambda[col.sel,,r]))/artefact | |
| 60 | print(max(delta)) | |
| 61 | densite = densite + piLambda[r] * | |
| 62 | det(rhoLambda[,,r])/(sqrt(2*base::pi))^m * exp(-tcrossprod(delta)/2.0) | |
| 63 | } | |
| 64 | llhLambda = c( sum(artefact^2 * log(densite)), (dimension+m+1)*k-1 ) | |
| 65 | list("phi"= phiLambda, "rho"= rhoLambda, "pi"= piLambda, "llh" = llhLambda) | |
| 66 | } | |
| 67 | ||
| 68 | # For each lambda, computation of the parameters | |
| 69 | out = | |
| 70 | if (ncores > 1) | |
| 71 | parLapply(cl, 1:length(S), computeAtLambda) | |
| 72 | else | |
| 73 | lapply(1:length(S), computeAtLambda) | |
| 74 | ||
| 75 | if (ncores > 1) | |
| 76 | parallel::stopCluster(cl) | |
| 77 | ||
| 78 | out | |
| c3bc4705 | 79 | } |