R/main.R

   1 SelMix = setRefClass(
   2     Class = "selmix",
   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 = "double",
  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 = "double",
  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,
  35         # initialisation for the reparametrized variance parameter
  36         rhoInit,
  37         # initialisation for the proportions
  38         piInit,
  39         # initialisation for the allocations probabilities in each component
  40         tauInit,
  41         # values for the regularization parameter grid
  42         gridLambda = c(),
  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,
  46         A2,
  47         # collection of estimations for the reparametrized conditional mean parameters
  48         Phi,
  49         # collection of estimations for the reparametrized variance parameters
  50         Rho,
  51         # collection of estimations for the proportions parameters
  52         Pi,
  53
  54         #immutable
  55         seuil = 1e-15
  56     ),
  57
  58     methods = list(
  59         #######################
  60         #initialize main object
  61         #######################
  62         initialize = function(X,Y,...)
  63         {
  64             "Initialize SelMix 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         },
  79
  80         ##################################
  81         #core workflow: compute all models
  82         ##################################
  83
  84         initParameters = function(k)
  85         {
  86             "Parameters initialization"
  87
  88             #smallEM initializes parameters by k-means and regression model in each component,
  89             #doing this 20 times, and keeping the values maximizing the likelihood after 10
  90             #iterations of the EM algorithm.
  91             init = initSmallEM(k,sx.X,sx.Y,sx.eps);
  92             phiInit <<- init$phi0;
  93             rhoInit <<- init$rho0;
  94             piInit  <<- init$pi0;
  95             tauInit <<- init$tau0;
  96         },
  97
  98         computeGridLambda = function()
  99         {
 100             "computation of the regularization grid"
 101             #(according to explicit formula given by EM algorithm)
 102
 103             gridLambda <<- gridLambda(sx.phiInit,sx.rhoInit,sx.piInit,sx.tauInit,sx.X,sx.Y,
 104                 sx.gamma,sx.mini,sx.maxi,sx.eps);
 105         },
 106
 107         computeRelevantParameters = function()
 108         {
 109             "Compute relevant parameters"
 110
 111             #select variables according to each regularization parameter
 112             #from the grid: sx.A1 corresponding to selected variables, and
 113             #sx.A2 corresponding to unselected variables.
 114             params = selectiontotale(sx.phiInit,sx.rhoInit,sx.piInit,sx.tauInit,sx.mini,sx.maxi,
 115                 sx.gamma,sx.gridLambda,sx.X,sx.Y,sx.seuil,sx.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             res = constructionModelesLassoMLE(sx.phiInit,sx.rhoInit,sx.piInit,sx.tauInit,
 129                 sx.mini,sx.maxi,sx.gamma,sx.gridLambda,sx.X,sx.Y,sx.seuil,sx.eps,sx.A1,sx.A2);
 130             return (list( phi=res$phi, rho=res$rho, pi=res$pi))
 131         },
 132
 133         runProcedure2 = function()
 134         {
 135             "Run procedure 2 [EMGrank]"
 136
 137             #compute parameter estimations, with the Low Rank
 138             #Estimator, restricted on selected variables.
 139             return (constructionModelesLassoRank(sx.Pi,sx.Rho,sx.mini,sx.maxi,sx.X,sx.Y,sx.eps,
 140                 sx.A1,sx.rangmin,sx.rangmax)$phi)
 141         },
 142
 143         run = function(procedure)
 144         {
 145             "main loop: over all k and all lambda"
 146
 147             # Run the all procedure, 1 with the
 148             #maximum likelihood refitting, and 2 with the Low Rank refitting.
 149                 p = dim(phiInit)[1]
 150                 m = dim(phiInit)[2]
 151                 for (k in kmin:kmax)
 152                 {
 153                     print(k)
 154                     initParameters(k)
 155                     computeGridLambda()
 156                     computeRelevantParameters()
 157                     if (procedure == 1)
 158                     {
 159                         r1 = runProcedure1(sx)
 160                         Phi2 = Phi
 161                         Rho2 = Rho
 162                         Pi2 = Pi
 163                         p = ncol(X)
 164                         m = ncol(Y)
 165                         if size(Phi2) == 0 #TODO: continue translation MATLAB --> R
 166                         Phi(:,:,1:k,:) = r1$phi;
 167                         Rho(:,:,1:k,:) = r1$rho;
 168                         Pi(1:k,:) = r1$pi;
 169                         else
 170                         Phi = zeros(p,m,sx.kmax,size(Phi2,4)+size(r1$phi,4));
 171                         Phi(:,:,1:size(Phi2,3),1:size(Phi2,4)) = Phi2;
 172                         Phi(:,:,1:k,size(Phi2,4)+1:end) = r1$phi;
 173                         Rho = zeros(m,m,sx.kmax,size(Rho2,4)+size(r1$rho,4));
 174                         Rho(:,:,1:size(Rho2,3),1:size(Rho2,4)) = Rho2;
 175                         Rho(:,:,1:k,size(Rho2,4)+1:end) = r1$rho;
 176                         Pi = zeros(sx.kmax,size(Pi2,2)+size(r1$pi,2));
 177                         Pi(1:size(Pi2,1),1:size(Pi2,2)) = Pi2;
 178                         Pi(1:k,size(Pi2,2)+1:end) = r1$pi;
 179                         end
 180                 else
 181                         [phi] = runProcedure2(sx);
 182                         phi
 183                         Phi2 = sx.Phi;
 184                         if size(Phi2,1) == 0
 185                         sx.Phi(:,:,1:k,:) = phi;
 186                         else
 187                         size(Phi2)
 188                         sx.Phi = zeros(p,m,sx.kmax,size(Phi2,4)+size(phi,4));
 189                         size(sx.Phi)
 190                         sx.Phi(:,:,1:size(Phi2,3),1:size(Phi2,4)) = Phi2;
 191                         sx.Phi(:,:,1:k,size(Phi2,4)+1:end) = phi;
 192                         end
 193
 194                 end
 195
 196
 197                 end
 198         end
 199
 200         ##################################################
 201         #pruning: select only one (or a few best ?!) model
 202         ##################################################
 203         #
 204         #       function[model] selectModel(sx)
 205         #           #TODO
 206         #           #model = sxModel(...);
 207         #       end
 208
 209         )
 210 )