fix typos
[valse.git] / R / main.R
CommitLineData
8e92c49c
BA
1Valse = setRefClass(
2 Class = "Valse",
39046da6
BA
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)
aa8df014 9 X = "matrix",
8e92c49c 10 # response data (size n*m, where n is the number of observations,
39046da6 11 # and m is the number of responses)
aa8df014 12 Y = "matrix",
39046da6
BA
13
14 # Optionally user defined (some default values)
15
16 # power in the penalty
09ab3c16 17 gamma = "numeric",
39046da6
BA
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
09ab3c16 23 eps = "numeric",
39046da6
BA
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
09ab3c16 34 phiInit = "numeric",
39046da6 35 # initialisation for the reparametrized variance parameter
09ab3c16 36 rhoInit = "numeric",
39046da6 37 # initialisation for the proportions
09ab3c16 38 piInit = "numeric",
39046da6 39 # initialisation for the allocations probabilities in each component
09ab3c16 40 tauInit = "numeric",
39046da6 41 # values for the regularization parameter grid
09ab3c16 42 gridLambda = "numeric",
39046da6
BA
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 ...
09ab3c16
BA
45 A1 = "integer",
46 A2 = "integer",
39046da6 47 # collection of estimations for the reparametrized conditional mean parameters
09ab3c16 48 Phi = "numeric",
39046da6 49 # collection of estimations for the reparametrized variance parameters
09ab3c16 50 Rho = "numeric",
39046da6 51 # collection of estimations for the proportions parameters
09ab3c16 52 Pi = "numeric",
39046da6 53
09ab3c16
BA
54 #immutable (TODO:?)
55 seuil = "numeric"
39046da6
BA
56 ),
57
58 methods = list(
59 #######################
60 #initialize main object
61 #######################
62 initialize = function(X,Y,...)
63 {
8e92c49c 64 "Initialize Valse object"
39046da6
BA
65
66 callSuper(...)
67
8e92c49c
BA
68 X <<- X
69 Y <<- Y
39046da6
BA
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))
09ab3c16 78 seuil <<- 1e-15 #immutable (TODO:?)
39046da6
BA
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.
8e92c49c
BA
92 init = initSmallEM(k,X,Y,eps)
93 phiInit <<- init$phi0
94 rhoInit <<- init$rho0
95 piInit <<- init$pi0
96 tauInit <<- init$tau0
39046da6
BA
97 },
98
99 computeGridLambda = function()
100 {
101 "computation of the regularization grid"
102 #(according to explicit formula given by EM algorithm)
103
8e92c49c 104 gridLambda <<- gridLambda(phiInit,rhoInit,piInit,tauInit,X,Y,gamma,mini,maxi,eps)
39046da6
BA
105 },
106
107 computeRelevantParameters = function()
108 {
109 "Compute relevant parameters"
110
111 #select variables according to each regularization parameter
8e92c49c
BA
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)
39046da6
BA
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.
8e92c49c
BA
128 return ( constructionModelesLassoMLE(
129 phiInit,rhoInit,piInit,tauInit,mini,maxi,gamma,gridLambda,X,Y,seuil,eps,A1,A2) )
39046da6
BA
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.
8e92c49c
BA
138 return ( constructionModelesLassoRank(Pi,Rho,mini,maxi,X,Y,eps,
139 A1,rangmin,rangmax) )
39046da6
BA
140 },
141
8e92c49c 142 run = function()
39046da6
BA
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.
8e92c49c
BA
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)
39046da6 157 {
8e92c49c
BA
158 r1 = runProcedure1()
159 Phi2 = Phi
160 Rho2 = Rho
161 Pi2 = Pi
162 p = ncol(X)
163 m = ncol(Y)
09ab3c16 164 if (is.null(dim(Phi2))) #test was: size(Phi2) == 0
8e92c49c 165 {
09ab3c16
BA
166 Phi[,,1:k] <<- r1$phi
167 Rho[,,1:k] <<- r1$rho
168 Pi[1:k,] <<- r1$pi
8e92c49c
BA
169 } else
170 {
09ab3c16
BA
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
8e92c49c
BA
180 }
181 } else
182 {
183 phi = runProcedure2()$phi
184 Phi2 = Phi
185 if (dim(Phi2)[1] == 0)
186 {
09ab3c16 187 Phi[,,1:k,] <<- phi
8e92c49c 188 } else
39046da6 189 {
09ab3c16
BA
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
8e92c49c
BA
193 }
194 }
195 }
196 }
197
39046da6 198 ##################################################
8e92c49c 199 #TODO: pruning: select only one (or a few best ?!) model
39046da6
BA
200 ##################################################
201 #
8e92c49c 202 # function[model] selectModel(
39046da6 203 # #TODO
8e92c49c 204 # #model = odel(...)
39046da6 205 # end
8e92c49c 206
39046da6
BA
207 )
208)