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