pkg/R/selectVariables.R

   1 #' selectVariables
   2 #' It is a function which construct, for a given lambda, the sets of relevant variables.
   3 #'
   4 #' @param phiInit an initial estimator for phi (size: p*m*k)
   5 #' @param rhoInit an initial estimator for rho (size: m*m*k)
   6 #' @param piInit    an initial estimator for pi (size : k)
   7 #' @param gamInit an initial estimator for gamma
   8 #' @param mini      minimum number of iterations in EM algorithm
   9 #' @param maxi      maximum number of iterations in EM algorithm
  10 #' @param gamma  power in the penalty
  11 #' @param glambda grid of regularization parameters
  12 #' @param X          matrix of regressors
  13 #' @param Y          matrix of responses
  14 #' @param thres  threshold to consider a coefficient to be equal to 0
  15 #' @param tau        threshold to say that EM algorithm has converged
  16 #'
  17 #' @return a list of outputs, for each lambda in grid: selected,Rho,Pi
  18 #'
  19 #' @examples TODO
  20 #'
  21 #' @export
  22 selectVariables = function(phiInit,rhoInit,piInit,gamInit,mini,maxi,gamma,glambda,X,Y,seuil,tau)
  23 {
  24     #TODO: parameter ncores (chaque tâche peut aussi demander du parallélisme...)
  25     cl = parallel::makeCluster( parallel::detectCores() / 4 )
  26     parallel::clusterExport(cl=cl,
  27         varlist=c("phiInit","rhoInit","gamInit","mini","maxi","glambda","X","Y","seuil","tau"),
  28         envir=environment())
  29     #Pour chaque lambda de la grille, on calcule les coefficients
  30     out = parLapply( seq_along(glambda), function(lambdaindex)
  31     {
  32         p = dim(phiInit)[1]
  33         m = dim(phiInit)[2]
  34
  35         params = EMGLLF(phiInit,rhoInit,piInit,gamInit,mini,maxi,gamma,glambda[lambdaIndex],X,Y,tau)
  36
  37         #selectedVariables: list where element j contains vector of selected variables in [1,m]
  38         selectedVariables = lapply(1:p, function(j) {
  39             #from boolean matrix mxk of selected variables obtain the corresponding boolean m-vector,
  40             #and finally return the corresponding indices
  41             seq_len(m)[ apply( abs(params$phi[j,,]) > seuil, 1, any ) ]
  42         })
  43
  44         list("selected"=selectedVariables,"Rho"=params$Rho,"Pi"=params$Pi)
  45     })
  46     parallel::stopCluster(cl)
  47     out
  48 }