R/main.R

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