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