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