X-Git-Url: https://git.auder.net/?p=valse.git;a=blobdiff_plain;f=pkg%2FR%2FconstructionModelesLassoMLE.R;h=1275ca3cc5d94c754ea3c741f605e78ac939e23a;hp=760da40be6952bad6312c36636d99df410009583;hb=ea5860f1b4fc91f06e371a0b26915198474a849d;hpb=ffdf94474d96cdd3e9d304ce809df7e62aa957ed diff --git a/pkg/R/constructionModelesLassoMLE.R b/pkg/R/constructionModelesLassoMLE.R index 760da40..1275ca3 100644 --- a/pkg/R/constructionModelesLassoMLE.R +++ b/pkg/R/constructionModelesLassoMLE.R @@ -21,8 +21,8 @@ #' #' @export constructionModelesLassoMLE <- function(phiInit, rhoInit, piInit, gamInit, mini, - maxi, gamma, X, Y, eps, S, ncores = 3, fast = TRUE, verbose = FALSE) - { + maxi, gamma, X, Y, eps, S, ncores = 3, fast, verbose) +{ if (ncores > 1) { cl <- parallel::makeCluster(ncores, outfile = "") @@ -30,62 +30,64 @@ constructionModelesLassoMLE <- function(phiInit, rhoInit, piInit, gamInit, mini, "rhoInit", "gamInit", "mini", "maxi", "gamma", "X", "Y", "eps", "S", "ncores", "fast", "verbose")) } - + # Individual model computation computeAtLambda <- function(lambda) { if (ncores > 1) require("valse") #nodes start with an empty environment - + if (verbose) print(paste("Computations for lambda=", lambda)) - - n <- dim(X)[1] - p <- dim(phiInit)[1] - m <- dim(phiInit)[2] - k <- dim(phiInit)[3] + + n <- nrow(X) + p <- ncol(X) + m <- ncol(Y) + k <- length(piInit) sel.lambda <- S[[lambda]]$selected # col.sel = which(colSums(sel.lambda)!=0) #if boolean matrix col.sel <- which(sapply(sel.lambda, length) > 0) #if list of selected vars if (length(col.sel) == 0) return(NULL) - + # lambda == 0 because we compute the EMV: no penalization here - res <- EMGLLF(phiInit[col.sel, , ], rhoInit, piInit, gamInit, mini, maxi, - gamma, 0, X[, col.sel], Y, eps, fast) - + res <- EMGLLF(array(phiInit,dim=c(p,m,k))[col.sel, , ], rhoInit, piInit, gamInit, + mini, maxi, gamma, 0, as.matrix(X[, col.sel]), Y, eps, fast) + # Eval dimension from the result + selected phiLambda2 <- res$phi rhoLambda <- res$rho piLambda <- res$pi phiLambda <- array(0, dim = c(p, m, k)) - for (j in seq_along(col.sel)) phiLambda[col.sel[j], sel.lambda[[j]], ] <- phiLambda2[j, - sel.lambda[[j]], ] + for (j in seq_along(col.sel)) + phiLambda[col.sel[j], sel.lambda[[j]], ] <- phiLambda2[j, sel.lambda[[j]], ] dimension <- length(unlist(sel.lambda)) - + # Computation of the loglikelihood densite <- vector("double", n) for (r in 1:k) { if (length(col.sel) == 1) { - delta <- (Y %*% rhoLambda[, , r] - (X[, col.sel] %*% t(phiLambda[col.sel, - , r]))) - } else delta <- (Y %*% rhoLambda[, , r] - (X[, col.sel] %*% phiLambda[col.sel, - , r])) - densite <- densite + piLambda[r] * det(rhoLambda[, , r])/(sqrt(2 * base::pi))^m * + delta <- (Y %*% rhoLambda[, , r] - (X[, col.sel] %*% t(phiLambda[col.sel, , r]))) + } else delta <- (Y %*% rhoLambda[, , r] - (X[, col.sel] %*% phiLambda[col.sel, , r])) + densite <- densite + piLambda[r] * gdet(rhoLambda[, , r])/(sqrt(2 * base::pi))^m * exp(-diag(tcrossprod(delta))/2) } llhLambda <- c(sum(log(densite)), (dimension + m + 1) * k - 1) list(phi = phiLambda, rho = rhoLambda, pi = piLambda, llh = llhLambda) } - + # For each lambda, computation of the parameters - out <- if (ncores > 1) - parLapply(cl, 1:length(S), computeAtLambda) else lapply(1:length(S), computeAtLambda) - + out <- + if (ncores > 1) { + parLapply(cl, 1:length(S), computeAtLambda) + } else { + lapply(1:length(S), computeAtLambda) + } + if (ncores > 1) parallel::stopCluster(cl) - + out }