Commit | Line | Data |
---|---|---|
4263503b | 1 | #' Wrapper function for OptimParams class |
cbd88fe5 | 2 | #' |
4263503b BA |
3 | #' @param X Data matrix of covariables |
4 | #' @param Y Output as a binary vector | |
2b3a6af5 BA |
5 | #' @param K Number of populations. |
6 | #' @param link The link type, 'logit' or 'probit'. | |
7 | #' @param M the empirical cross-moments between X and Y (optional) | |
cbd88fe5 | 8 | #' |
2b3a6af5 BA |
9 | #' @return An object 'op' of class OptimParams, initialized so that |
10 | #' \code{op$run(θ0)} outputs the list of optimized parameters | |
cbd88fe5 BA |
11 | #' \itemize{ |
12 | #' \item p: proportions, size K | |
13 | #' \item β: regression matrix, size dxK | |
14 | #' \item b: intercepts, size K | |
15 | #' } | |
2b3a6af5 BA |
16 | #' θ0 is a list containing the initial parameters. Only β is required |
17 | #' (p would be set to (1/K,...,1/K) and b to (0,...0)). | |
cbd88fe5 BA |
18 | #' |
19 | #' @seealso \code{multiRun} to estimate statistics based on β, and | |
20 | #' \code{generateSampleIO} for I/O random generation. | |
21 | #' | |
22 | #' @examples | |
23 | #' # Optimize parameters from estimated μ | |
2b3a6af5 BA |
24 | #' io <- generateSampleIO(100, |
25 | #' 1/2, matrix(c(1,-2,3,1),ncol=2), c(0,0), "logit") | |
cbd88fe5 | 26 | #' μ = computeMu(io$X, io$Y, list(K=2)) |
4263503b | 27 | #' o <- optimParams(io$X, io$Y, 2, "logit") |
2b3a6af5 | 28 | #' \donttest{ |
7737c2fa BA |
29 | #' θ0 <- list(p=1/2, β=μ, b=c(0,0)) |
30 | #' par0 <- o$run(θ0) | |
cbd88fe5 | 31 | #' # Compare with another starting point |
7737c2fa BA |
32 | #' θ1 <- list(p=1/2, β=2*μ, b=c(0,0)) |
33 | #' par1 <- o$run(θ1) | |
2b3a6af5 | 34 | #' # Look at the function values at par0 and par1: |
cbd88fe5 | 35 | #' o$f( o$linArgs(par0) ) |
2b3a6af5 BA |
36 | #' o$f( o$linArgs(par1) )} |
37 | #' | |
cbd88fe5 | 38 | #' @export |
f4e42a2b | 39 | optimParams <- function(X, Y, K, link=c("logit","probit"), M=NULL) |
cbd88fe5 | 40 | { |
6dd5c2ac | 41 | # Check arguments |
4263503b BA |
42 | if (!is.matrix(X) || any(is.na(X))) |
43 | stop("X: numeric matrix, no NAs") | |
0a630686 | 44 | if (!is.numeric(Y) || any(is.na(Y)) || any(Y!=0 & Y!=1)) |
4263503b | 45 | stop("Y: binary vector with 0 and 1 only") |
6dd5c2ac | 46 | link <- match.arg(link) |
4263503b BA |
47 | if (!is.numeric(K) || K!=floor(K) || K < 2) |
48 | stop("K: integer >= 2") | |
cbd88fe5 | 49 | |
f4e42a2b BA |
50 | if (is.null(M)) |
51 | { | |
52 | # Precompute empirical moments | |
53 | Mtmp <- computeMoments(X, Y) | |
54 | M1 <- as.double(Mtmp[[1]]) | |
55 | M2 <- as.double(Mtmp[[2]]) | |
56 | M3 <- as.double(Mtmp[[3]]) | |
57 | M <- c(M1, M2, M3) | |
58 | } | |
d9edcd6c BA |
59 | else |
60 | M <- c(M[[1]], M[[2]], M[[3]]) | |
f4e42a2b | 61 | |
6dd5c2ac BA |
62 | # Build and return optimization algorithm object |
63 | methods::new("OptimParams", "li"=link, "X"=X, | |
f4e42a2b | 64 | "Y"=as.integer(Y), "K"=as.integer(K), "Mhat"=as.double(M)) |
cbd88fe5 BA |
65 | } |
66 | ||
2b3a6af5 BA |
67 | # Encapsulated optimization for p (proportions), β and b (regression parameters) |
68 | # | |
69 | # Optimize the parameters of a mixture of logistic regressions model, possibly using | |
70 | # \code{mu <- computeMu(...)} as a partial starting point. | |
71 | # | |
72 | # @field li Link function, 'logit' or 'probit' | |
73 | # @field X Data matrix of covariables | |
74 | # @field Y Output as a binary vector | |
75 | # @field Mhat Vector of empirical moments | |
76 | # @field K Number of populations | |
77 | # @field n Number of sample points | |
78 | # @field d Number of dimensions | |
79 | # @field W Weights matrix (initialized at identity) | |
80 | # | |
cbd88fe5 | 81 | setRefClass( |
6dd5c2ac | 82 | Class = "OptimParams", |
cbd88fe5 | 83 | |
6dd5c2ac BA |
84 | fields = list( |
85 | # Inputs | |
86 | li = "character", #link function | |
87 | X = "matrix", | |
88 | Y = "numeric", | |
7737c2fa | 89 | Mhat = "numeric", #vector of empirical moments |
6dd5c2ac BA |
90 | # Dimensions |
91 | K = "integer", | |
4263503b | 92 | n = "integer", |
6dd5c2ac | 93 | d = "integer", |
e92d9d9d BA |
94 | # Weights matrix (generalized least square) |
95 | W = "matrix" | |
6dd5c2ac | 96 | ), |
cbd88fe5 | 97 | |
6dd5c2ac BA |
98 | methods = list( |
99 | initialize = function(...) | |
100 | { | |
101 | "Check args and initialize K, d, W" | |
cbd88fe5 | 102 | |
4263503b | 103 | callSuper(...) |
f4e42a2b BA |
104 | if (!hasArg("X") || !hasArg("Y") || !hasArg("K") |
105 | || !hasArg("li") || !hasArg("Mhat")) | |
106 | { | |
6dd5c2ac | 107 | stop("Missing arguments") |
f4e42a2b | 108 | } |
4263503b | 109 | |
6dd5c2ac | 110 | n <<- nrow(X) |
f4e42a2b | 111 | d <<- ncol(X) |
86f257f8 | 112 | # W will be initialized when calling run() |
6dd5c2ac | 113 | }, |
cbd88fe5 | 114 | |
6dd5c2ac BA |
115 | expArgs = function(v) |
116 | { | |
117 | "Expand individual arguments from vector v into a list" | |
cbd88fe5 | 118 | |
6dd5c2ac BA |
119 | list( |
120 | # p: dimension K-1, need to be completed | |
121 | "p" = c(v[1:(K-1)], 1-sum(v[1:(K-1)])), | |
44559add | 122 | "β" = t(matrix(v[K:(K+d*K-1)], ncol=d)), |
6dd5c2ac BA |
123 | "b" = v[(K+d*K):(K+(d+1)*K-1)]) |
124 | }, | |
cbd88fe5 | 125 | |
6dd5c2ac BA |
126 | linArgs = function(L) |
127 | { | |
128 | "Linearize vectors+matrices from list L into a vector" | |
cbd88fe5 | 129 | |
44559add BA |
130 | # β linearized row by row, to match derivatives order |
131 | c(L$p[1:(K-1)], as.double(t(L$β)), L$b) | |
6dd5c2ac | 132 | }, |
cbd88fe5 | 133 | |
2b3a6af5 | 134 | # TODO: relocate computeW in utils.R |
7737c2fa | 135 | computeW = function(θ) |
4263503b | 136 | { |
2b3a6af5 BA |
137 | "Compute the weights matrix from a parameters list" |
138 | ||
0f5fbd13 | 139 | require(MASS) |
4bf8494d | 140 | dd <- d + d^2 + d^3 |
9a6881ed BA |
141 | M <- Moments(θ) |
142 | Omega <- matrix( .C("Compute_Omega", | |
074c721a | 143 | X=as.double(X), Y=as.integer(Y), M=as.double(M), |
9a6881ed BA |
144 | pn=as.integer(n), pd=as.integer(d), |
145 | W=as.double(W), PACKAGE="morpheus")$W, nrow=dd, ncol=dd ) | |
0f5fbd13 | 146 | MASS::ginv(Omega) |
4263503b BA |
147 | }, |
148 | ||
b389a46a | 149 | Moments = function(θ) |
4263503b | 150 | { |
2b3a6af5 | 151 | "Compute the vector of theoretical moments (size d+d^2+d^3)" |
cbd88fe5 | 152 | |
7737c2fa | 153 | p <- θ$p |
6dd5c2ac BA |
154 | β <- θ$β |
155 | λ <- sqrt(colSums(β^2)) | |
156 | b <- θ$b | |
157 | ||
158 | # Tensorial products β^2 = β2 and β^3 = β3 must be computed from current β1 | |
159 | β2 <- apply(β, 2, function(col) col %o% col) | |
160 | β3 <- apply(β, 2, function(col) col %o% col %o% col) | |
161 | ||
162 | c( | |
163 | β %*% (p * .G(li,1,λ,b)), | |
164 | β2 %*% (p * .G(li,2,λ,b)), | |
165 | β3 %*% (p * .G(li,3,λ,b))) | |
7737c2fa BA |
166 | }, |
167 | ||
168 | f = function(θ) | |
169 | { | |
2b3a6af5 | 170 | "Function to minimize: t(hat_Mi - Mi(θ)) . W . (hat_Mi - Mi(θ))" |
7737c2fa | 171 | |
0a630686 | 172 | L <- expArgs(θ) |
6dd5c2ac | 173 | A <- as.matrix(Mhat - Moments(L)) |
4263503b BA |
174 | t(A) %*% W %*% A |
175 | }, | |
cbd88fe5 | 176 | |
6dd5c2ac BA |
177 | grad_f = function(θ) |
178 | { | |
2b3a6af5 | 179 | "Gradient of f: vector of size (K-1) + d*K + K = (d+2)*K - 1" |
cbd88fe5 | 180 | |
0a630686 | 181 | L <- expArgs(θ) |
0f5fbd13 | 182 | -2 * t(grad_M(L)) %*% W %*% as.matrix(Mhat - Moments(L)) |
b389a46a | 183 | }, |
4263503b | 184 | |
7737c2fa | 185 | grad_M = function(θ) |
4263503b | 186 | { |
2b3a6af5 | 187 | "Gradient of the moments vector: matrix of size d+d^2+d^3 x K-1+K+d*K" |
4263503b | 188 | |
6dd5c2ac BA |
189 | p <- θ$p |
190 | β <- θ$β | |
191 | λ <- sqrt(colSums(β^2)) | |
192 | μ <- sweep(β, 2, λ, '/') | |
193 | b <- θ$b | |
7737c2fa BA |
194 | |
195 | res <- matrix(nrow=nrow(W), ncol=0) | |
cbd88fe5 | 196 | |
6dd5c2ac BA |
197 | # Tensorial products β^2 = β2 and β^3 = β3 must be computed from current β1 |
198 | β2 <- apply(β, 2, function(col) col %o% col) | |
199 | β3 <- apply(β, 2, function(col) col %o% col %o% col) | |
cbd88fe5 | 200 | |
6dd5c2ac BA |
201 | # Some precomputations |
202 | G1 = .G(li,1,λ,b) | |
203 | G2 = .G(li,2,λ,b) | |
204 | G3 = .G(li,3,λ,b) | |
205 | G4 = .G(li,4,λ,b) | |
206 | G5 = .G(li,5,λ,b) | |
cbd88fe5 | 207 | |
7737c2fa | 208 | # Gradient on p: K-1 columns, dim rows |
6dd5c2ac BA |
209 | km1 = 1:(K-1) |
210 | res <- cbind(res, rbind( | |
0a630686 BA |
211 | sweep(as.matrix(β [,km1]), 2, G1[km1], '*') - G1[K] * β [,K], |
212 | sweep(as.matrix(β2[,km1]), 2, G2[km1], '*') - G2[K] * β2[,K], | |
213 | sweep(as.matrix(β3[,km1]), 2, G3[km1], '*') - G3[K] * β3[,K] )) | |
cbd88fe5 | 214 | |
6dd5c2ac BA |
215 | for (i in 1:d) |
216 | { | |
217 | # i determines the derivated matrix dβ[2,3] | |
218 | ||
219 | dβ_left <- sweep(β, 2, p * G3 * β[i,], '*') | |
220 | dβ_right <- matrix(0, nrow=d, ncol=K) | |
221 | block <- i | |
222 | dβ_right[block,] <- dβ_right[block,] + 1 | |
223 | dβ <- dβ_left + sweep(dβ_right, 2, p * G1, '*') | |
224 | ||
225 | dβ2_left <- sweep(β2, 2, p * G4 * β[i,], '*') | |
226 | dβ2_right <- do.call( rbind, lapply(1:d, function(j) { | |
227 | sweep(dβ_right, 2, β[j,], '*') | |
228 | }) ) | |
229 | block <- ((i-1)*d+1):(i*d) | |
230 | dβ2_right[block,] <- dβ2_right[block,] + β | |
231 | dβ2 <- dβ2_left + sweep(dβ2_right, 2, p * G2, '*') | |
232 | ||
233 | dβ3_left <- sweep(β3, 2, p * G5 * β[i,], '*') | |
234 | dβ3_right <- do.call( rbind, lapply(1:d, function(j) { | |
235 | sweep(dβ2_right, 2, β[j,], '*') | |
236 | }) ) | |
237 | block <- ((i-1)*d*d+1):(i*d*d) | |
238 | dβ3_right[block,] <- dβ3_right[block,] + β2 | |
239 | dβ3 <- dβ3_left + sweep(dβ3_right, 2, p * G3, '*') | |
240 | ||
241 | res <- cbind(res, rbind(dβ, dβ2, dβ3)) | |
242 | } | |
cbd88fe5 | 243 | |
7737c2fa | 244 | # Gradient on b |
6dd5c2ac BA |
245 | res <- cbind(res, rbind( |
246 | sweep(β, 2, p * G2, '*'), | |
247 | sweep(β2, 2, p * G3, '*'), | |
248 | sweep(β3, 2, p * G4, '*') )) | |
cbd88fe5 | 249 | |
6dd5c2ac BA |
250 | res |
251 | }, | |
cbd88fe5 | 252 | |
6dd5c2ac BA |
253 | run = function(θ0) |
254 | { | |
255 | "Run optimization from θ0 with solver..." | |
7737c2fa | 256 | |
6dd5c2ac BA |
257 | if (!is.list(θ0)) |
258 | stop("θ0: list") | |
7737c2fa BA |
259 | if (is.null(θ0$β)) |
260 | stop("At least θ0$β must be provided") | |
0f5fbd13 BA |
261 | if (!is.matrix(θ0$β) || any(is.na(θ0$β)) |
262 | || nrow(θ0$β) != d || ncol(θ0$β) != K) | |
263 | { | |
264 | stop("θ0$β: matrix, no NA, nrow = d, ncol = K") | |
265 | } | |
7737c2fa BA |
266 | if (is.null(θ0$p)) |
267 | θ0$p = rep(1/K, K-1) | |
0f5fbd13 BA |
268 | else if (!is.numeric(θ0$p) || length(θ0$p) != K-1 |
269 | || any(is.na(θ0$p)) || sum(θ0$p) > 1) | |
270 | { | |
271 | stop("θ0$p: length K-1, no NA, positive integers, sum to <= 1") | |
272 | } | |
273 | if (is.null(θ0$b)) | |
7737c2fa | 274 | θ0$b = rep(0, K) |
0f5fbd13 BA |
275 | else if (!is.numeric(θ0$b) || length(θ0$b) != K || any(is.na(θ0$b))) |
276 | stop("θ0$b: length K, no NA") | |
86f257f8 BA |
277 | |
278 | # (Re)Set W to identity, to allow several run from the same object | |
279 | W <<- diag(d+d^2+d^3) | |
280 | ||
07f2d045 | 281 | loopMax <- 2 #TODO: loopMax = 3 ? Seems not improving... |
2591fa83 | 282 | x_init <- linArgs(θ0) |
ef0d907c | 283 | for (loop in 1:loopMax) |
4bf8494d | 284 | { |
2591fa83 | 285 | op_res = constrOptim( x_init, .self$f, .self$grad_f, |
4bf8494d BA |
286 | ui=cbind( |
287 | rbind( rep(-1,K-1), diag(K-1) ), | |
288 | matrix(0, nrow=K, ncol=(d+1)*K) ), | |
289 | ci=c(-1,rep(0,K-1)) ) | |
ef0d907c BA |
290 | if (loop < loopMax) #avoid computing an extra W |
291 | W <<- computeW(expArgs(op_res$par)) | |
2989133a | 292 | #x_init <- op_res$par #degrades performances (TODO: why?) |
4bf8494d | 293 | } |
4263503b | 294 | |
6dd5c2ac BA |
295 | expArgs(op_res$par) |
296 | } | |
297 | ) | |
cbd88fe5 BA |
298 | ) |
299 | ||
300 | # Compute vectorial E[g^{(order)}(<β,x> + b)] with x~N(0,Id) (integral in R^d) | |
301 | # = E[g^{(order)}(z)] with z~N(b,diag(λ)) | |
4263503b | 302 | # by numerically evaluating the integral. |
cbd88fe5 BA |
303 | # |
304 | # @param link Link, 'logit' or 'probit' | |
305 | # @param order Order of derivative | |
306 | # @param λ Norm of columns of β | |
307 | # @param b Intercept | |
308 | # | |
309 | .G <- function(link, order, λ, b) | |
310 | { | |
6dd5c2ac BA |
311 | # NOTE: weird "integral divergent" error on inputs: |
312 | # link="probit"; order=2; λ=c(531.8099,586.8893,523.5816); b=c(-118.512674,-3.488020,2.109969) | |
313 | # Switch to pracma package for that (but it seems slow...) | |
4263503b BA |
314 | sapply( seq_along(λ), function(k) { |
315 | res <- NULL | |
316 | tryCatch({ | |
317 | # Fast code, may fail: | |
318 | res <- stats::integrate( | |
319 | function(z) .deriv[[link]][[order]](λ[k]*z+b[k]) * exp(-z^2/2) / sqrt(2*pi), | |
320 | lower=-Inf, upper=Inf )$value | |
321 | }, error = function(e) { | |
322 | # Robust slow code, no fails observed: | |
323 | sink("/dev/null") #pracma package has some useless printed outputs... | |
324 | res <- pracma::integral( | |
325 | function(z) .deriv[[link]][[order]](λ[k]*z+b[k]) * exp(-z^2/2) / sqrt(2*pi), | |
326 | xmin=-Inf, xmax=Inf, method="Kronrod") | |
327 | sink() | |
328 | }) | |
329 | res | |
330 | }) | |
cbd88fe5 BA |
331 | } |
332 | ||
333 | # Derivatives list: g^(k)(x) for links 'logit' and 'probit' | |
334 | # | |
335 | .deriv <- list( | |
6dd5c2ac BA |
336 | "probit"=list( |
337 | # 'probit' derivatives list; | |
338 | # NOTE: exact values for the integral E[g^(k)(λz+b)] could be computed | |
339 | function(x) exp(-x^2/2)/(sqrt(2*pi)), #g' | |
340 | function(x) exp(-x^2/2)/(sqrt(2*pi)) * -x, #g'' | |
341 | function(x) exp(-x^2/2)/(sqrt(2*pi)) * ( x^2 - 1), #g^(3) | |
342 | function(x) exp(-x^2/2)/(sqrt(2*pi)) * (-x^3 + 3*x), #g^(4) | |
343 | function(x) exp(-x^2/2)/(sqrt(2*pi)) * ( x^4 - 6*x^2 + 3) #g^(5) | |
344 | ), | |
345 | "logit"=list( | |
346 | # Sigmoid derivatives list, obtained with http://www.derivative-calculator.net/ | |
347 | # @seealso http://www.ece.uc.edu/~aminai/papers/minai_sigmoids_NN93.pdf | |
348 | function(x) {e=exp(x); .zin(e /(e+1)^2)}, #g' | |
349 | function(x) {e=exp(x); .zin(e*(-e + 1) /(e+1)^3)}, #g'' | |
350 | function(x) {e=exp(x); .zin(e*( e^2 - 4*e + 1) /(e+1)^4)}, #g^(3) | |
351 | function(x) {e=exp(x); .zin(e*(-e^3 + 11*e^2 - 11*e + 1) /(e+1)^5)}, #g^(4) | |
352 | function(x) {e=exp(x); .zin(e*( e^4 - 26*e^3 + 66*e^2 - 26*e + 1)/(e+1)^6)} #g^(5) | |
353 | ) | |
cbd88fe5 BA |
354 | ) |
355 | ||
356 | # Utility for integration: "[return] zero if [argument is] NaN" (Inf / Inf divs) | |
357 | # | |
358 | # @param x Ratio of polynoms of exponentials, as in .S[[i]] | |
359 | # | |
360 | .zin <- function(x) | |
361 | { | |
6dd5c2ac BA |
362 | x[is.nan(x)] <- 0. |
363 | x | |
cbd88fe5 | 364 | } |