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