| 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 = "matrix", |
| 10 | # response data (size n*m, where n is the number of observations, |
| 11 | # and m is the number of responses) |
| 12 | Y = "matrix", |
| 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 | ) |