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