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