9 # regression data (size n*p, where n is the number of observations,
10 # and p is the number of regressors)
12 # response data (size n*m, where n is the number of observations,
13 # and m is the number of responses)
16 # Optionally user defined (some default values)
18 # power in the penalty
20 # minimum number of iterations for EM algorithm
22 # maximum number of iterations for EM algorithm
24 # threshold for stopping EM algorithm
26 # minimum number of components in the mixture
28 # maximum number of components in the mixture
33 # Computed through the workflow
35 # initialisation for the reparametrized conditional mean parameter
37 # initialisation for the reparametrized variance parameter
39 # initialisation for the proportions
41 # initialisation for the allocations probabilities in each component
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 ...
49 # collection of estimations for the reparametrized conditional mean parameters
51 # collection of estimations for the reparametrized variance parameters
53 # collection of estimations for the proportions parameters
61 #######################
62 #initialize main object
63 #######################
64 initialize = function(X,Y,...)
66 "Initialize Valse object"
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:?)
83 ##################################
84 #core workflow: compute all models
85 ##################################
87 initParameters = function(k)
89 "Parameters initialization"
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)
101 computeGridLambda = function()
103 "computation of the regularization grid"
104 #(according to explicit formula given by EM algorithm)
106 gridLambda <<- gridLambda(phiInit,rhoInit,piInit,tauInit,X,Y,gamma,mini,maxi,eps)
109 computeRelevantParameters = function()
111 "Compute relevant parameters"
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)
124 runProcedure1 = function()
126 "Run procedure 1 [EMGLLF]"
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) )
134 runProcedure2 = function()
136 "Run procedure 2 [EMGrank]"
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) )
146 "main loop: over all k and all lambda"
148 # Run the all procedure, 1 with the
149 #maximum likelihood refitting, and 2 with the Low Rank refitting.
157 computeRelevantParameters()
166 if (is.null(dim(Phi2))) #test was: size(Phi2) == 0
168 Phi[,,1:k] <<- r1$phi
169 Rho[,,1:k] <<- r1$rho
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
185 phi = runProcedure2()$phi
187 if (dim(Phi2)[1] == 0)
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
200 ##################################################
201 #TODO: pruning: select only one (or a few best ?!) model
202 ##################################################
204 # function[model] selectModel(