pkg/R/selectVariables.R

   1 #' selectVariables
   2 #'
   3 #' It is a function which construct, for a given lambda, the sets of relevant variables.
   4 #'
   5 #' @param phiInit an initial estimator for phi (size: p*m*k)
   6 #' @param rhoInit an initial estimator for rho (size: m*m*k)
   7 #' @param piInit    an initial estimator for pi (size : k)
   8 #' @param gamInit an initial estimator for gamma
   9 #' @param mini      minimum number of iterations in EM algorithm
  10 #' @param maxi      maximum number of iterations in EM algorithm
  11 #' @param gamma  power in the penalty
  12 #' @param glambda grid of regularization parameters
  13 #' @param X          matrix of regressors
  14 #' @param Y          matrix of responses
  15 #' @param thres  threshold to consider a coefficient to be equal to 0
  16 #' @param tau        threshold to say that EM algorithm has converged
  17 #' @param ncores Number or cores for parallel execution (1 to disable)
  18 #'
  19 #' @return a list of outputs, for each lambda in grid: selected,Rho,Pi
  20 #'
  21 #' @examples TODO
  22 #'
  23 #' @export
  24 #'
  25 selectVariables = function(phiInit,rhoInit,piInit,gamInit,mini,maxi,gamma,glambda,
  26     X,Y,thresh,tau, ncores=3)
  27 {
  28     if (ncores > 1)
  29     {
  30         cl = parallel::makeCluster(ncores)
  31         parallel::clusterExport(cl=cl,
  32             varlist=c("phiInit","rhoInit","gamInit","mini","maxi","glambda","X","Y","thresh","tau"),
  33             envir=environment())
  34     }
  35
  36     # Calcul pour un lambda
  37     computeCoefs <-function(lambda)
  38     {
  39         params = EMGLLF(phiInit,rhoInit,piInit,gamInit,mini,maxi,gamma,lambda,X,Y,tau)
  40
  41         p = dim(phiInit)[1]
  42         m = dim(phiInit)[2]
  43
  44         #selectedVariables: list where element j contains vector of selected variables in [1,m]
  45         selectedVariables = lapply(1:p, function(j) {
  46             #from boolean matrix mxk of selected variables obtain the corresponding boolean m-vector,
  47             #and finally return the corresponding indices
  48           seq_len(m)[ apply( abs(params$phi[j,,]) > thresh, 1, any ) ]
  49         })
  50
  51         list("selected"=selectedVariables,"Rho"=params$rho,"Pi"=params$pi)
  52     }
  53
  54     # Pour chaque lambda de la grille, on calcule les coefficients
  55     out <-
  56         if (ncores > 1)
  57             parLapply(cl, glambda, computeCoefs)
  58         else
  59             lapply(glambda, computeCoefs)
  60     if (ncores > 1)
  61         parallel::stopCluster(cl)
  62
  63     # Suppression doublons
  64     sha1_array <- lapply(out, digest::sha1)
  65     out[ !duplicated(sha1_array) ]
  66
  67     out
  68 }