R/main.R

   1 #' @useDynLib valse
   2
   3 Valse = setRefClass(
   4     Class = "Valse",
   5
   6     fields = c(
   7         # User defined
   8
   9         # regression data (size n*p, where n is the number of observations,
  10         # and p is the number of regressors)
  11         X = "matrix",
  12         # response data (size n*m, where n is the number of observations,
  13         # and m is the number of responses)
  14         Y = "matrix",
  15
  16         # Optionally user defined (some default values)
  17
  18         # power in the penalty
  19         gamma = "numeric",
  20         # minimum number of iterations for EM algorithm
  21         mini = "integer",
  22         # maximum number of iterations for EM algorithm
  23         maxi = "integer",
  24         # threshold for stopping EM algorithm
  25         eps = "numeric",
  26         # minimum number of components in the mixture
  27         kmin = "integer",
  28         # maximum number of components in the mixture
  29         kmax = "integer",
  30         rangmin = "integer",
  31         rangmax = "integer",
  32
  33         # Computed through the workflow
  34
  35         # initialisation for the reparametrized conditional mean parameter
  36         phiInit = "numeric",
  37         # initialisation for the reparametrized variance parameter
  38         rhoInit = "numeric",
  39         # initialisation for the proportions
  40         piInit = "numeric",
  41         # initialisation for the allocations probabilities in each component
  42         tauInit = "numeric",
  43         # values for the regularization parameter grid
  44         gridLambda = "numeric",
  45         # je ne crois pas vraiment qu'il faille les mettre en sortie, d'autant plus qu'on construit
  46         # une matrice A1 et A2 pour chaque k, et elles sont grandes, donc ca coute un peu cher ...
  47         A1 = "integer",
  48         A2 = "integer",
  49         # collection of estimations for the reparametrized conditional mean parameters
  50         Phi = "numeric",
  51         # collection of estimations for the reparametrized variance parameters
  52         Rho = "numeric",
  53         # collection of estimations for the proportions parameters
  54         Pi = "numeric",
  55
  56         #immutable (TODO:?)
  57         seuil = "numeric"
  58     ),
  59
  60     methods = list(
  61         #######################
  62         #initialize main object
  63         #######################
  64         initialize = function(X,Y,...)
  65         {
  66             "Initialize Valse object"
  67
  68             callSuper(...)
  69
  70             X <<- X
  71             Y <<- Y
  72             gamma <<- ifelse (hasArg("gamma"), gamma, 1.)
  73             mini <<- ifelse (hasArg("mini"), mini, as.integer(5))
  74             maxi <<- ifelse (hasArg("maxi"), maxi, as.integer(10))
  75             eps <<- ifelse (hasArg("eps"), eps, 1e-6)
  76             kmin <<- ifelse (hasArg("kmin"), kmin, as.integer(2))
  77             kmax <<- ifelse (hasArg("kmax"), kmax, as.integer(3))
  78             rangmin <<- ifelse (hasArg("rangmin"), rangmin, as.integer(2))
  79             rangmax <<- ifelse (hasArg("rangmax"), rangmax, as.integer(3))
  80             seuil <<- 1e-15 #immutable (TODO:?)
  81         },
  82
  83         ##################################
  84         #core workflow: compute all models
  85         ##################################
  86
  87         initParameters = function(k)
  88         {
  89             "Parameters initialization"
  90
  91             #smallEM initializes parameters by k-means and regression model in each component,
  92             #doing this 20 times, and keeping the values maximizing the likelihood after 10
  93             #iterations of the EM algorithm.
  94             init = initSmallEM(k,X,Y,eps)
  95             phiInit <<- init$phi0
  96             rhoInit <<- init$rho0
  97             piInit  <<- init$pi0
  98             tauInit <<- init$tau0
  99         },
 100
 101         computeGridLambda = function()
 102         {
 103             "computation of the regularization grid"
 104             #(according to explicit formula given by EM algorithm)
 105
 106             gridLambda <<- gridLambda(phiInit,rhoInit,piInit,tauInit,X,Y,gamma,mini,maxi,eps)
 107         },
 108
 109         computeRelevantParameters = function()
 110         {
 111             "Compute relevant parameters"
 112
 113             #select variables according to each regularization parameter
 114             #from the grid: A1 corresponding to selected variables, and
 115             #A2 corresponding to unselected variables.
 116             params = selectiontotale(
 117                 phiInit,rhoInit,piInit,tauInit,mini,maxi,gamma,gridLambda,X,Y,seuil,eps)
 118             A1 <<- params$A1
 119             A2 <<- params$A2
 120             Rho <<- params$Rho
 121             Pi <<- params$Pi
 122         },
 123
 124         runProcedure1 = function()
 125         {
 126             "Run procedure 1 [EMGLLF]"
 127
 128             #compute parameter estimations, with the Maximum Likelihood
 129             #Estimator, restricted on selected variables.
 130             return ( constructionModelesLassoMLE(
 131                 phiInit,rhoInit,piInit,tauInit,mini,maxi,gamma,gridLambda,X,Y,seuil,eps,A1,A2) )
 132         },
 133
 134         runProcedure2 = function()
 135         {
 136             "Run procedure 2 [EMGrank]"
 137
 138             #compute parameter estimations, with the Low Rank
 139             #Estimator, restricted on selected variables.
 140             return ( constructionModelesLassoRank(Pi,Rho,mini,maxi,X,Y,eps,
 141                 A1,rangmin,rangmax) )
 142         },
 143
 144         run = function()
 145         {
 146             "main loop: over all k and all lambda"
 147
 148             # Run the all procedure, 1 with the
 149             #maximum likelihood refitting, and 2 with the Low Rank refitting.
 150             p = dim(phiInit)[1]
 151             m = dim(phiInit)[2]
 152             for (k in kmin:kmax)
 153             {
 154                 print(k)
 155                 initParameters(k)
 156                 computeGridLambda()
 157                 computeRelevantParameters()
 158                 if (procedure == 1)
 159                 {
 160                     r1 = runProcedure1()
 161                     Phi2 = Phi
 162                     Rho2 = Rho
 163                     Pi2 = Pi
 164                     p = ncol(X)
 165                     m = ncol(Y)
 166                     if (is.null(dim(Phi2))) #test was: size(Phi2) == 0
 167                     {
 168                         Phi[,,1:k] <<- r1$phi
 169                         Rho[,,1:k] <<- r1$rho
 170                         Pi[1:k,] <<- r1$pi
 171                     } else
 172                     {
 173                         Phi <<- array(0., dim=c(p,m,kmax,dim(Phi2)[4]+dim(r1$phi)[4]))
 174                         Phi[,,1:(dim(Phi2)[3]),1:(dim(Phi2)[4])] <<- Phi2
 175                         Phi[,,1:k,dim(Phi2)[4]+1] <<- r1$phi
 176                         Rho <<- array(0., dim=c(m,m,kmax,dim(Rho2)[4]+dim(r1$rho)[4]))
 177                         Rho[,,1:(dim(Rho2)[3]),1:(dim(Rho2)[4])] <<- Rho2
 178                         Rho[,,1:k,dim(Rho2)[4]+1] <<- r1$rho
 179                         Pi <<- array(0., dim=c(kmax,dim(Pi2)[2]+dim(r1$pi)[2]))
 180                         Pi[1:nrow(Pi2),1:ncol(Pi2)] <<- Pi2
 181                         Pi[1:k,ncol(Pi2)+1] <<- r1$pi
 182                     }
 183                 } else
 184                 {
 185                     phi = runProcedure2()$phi
 186                     Phi2 = Phi
 187                     if (dim(Phi2)[1] == 0)
 188                     {
 189                         Phi[,,1:k,] <<- phi
 190                     } else
 191                     {
 192                         Phi <<- array(0., dim=c(p,m,kmax,dim(Phi2)[4]+dim(phi)[4]))
 193                         Phi[,,1:(dim(Phi2)[3]),1:(dim(Phi2)[4])] <<- Phi2
 194                         Phi[,,1:k,-(1:(dim(Phi2)[4]))] <<- phi
 195                     }
 196                 }
 197             }
 198         }
 199
 200         ##################################################
 201         #TODO: pruning: select only one (or a few best ?!) model
 202         ##################################################
 203         #
 204         #       function[model] selectModel(
 205         #           #TODO
 206         #           #model = odel(...)
 207         #       end
 208
 209         )
 210 )