Fix Y computation in data sampling
[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 #' \item weights Weights on moments when minimizing sum of squares
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
60 weights <- optargs$weights
61 if (is.null(weights))
62 weights <- rep(1, K)
63
64 # Build and return optimization algorithm object
65 methods::new("OptimParams", "li"=link, "M1"=as.double(M[[1]]),
66 "M2"=as.double(M[[2]]), "M3"=as.double(M[[3]]),
67 "weights"=weights, "K"=as.integer(K))
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)
76 # @field weights Vector of moments' weights
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
89 weights = "numeric", #weights on moments
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(
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 ) )
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
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
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
208 grad_β[i,] <-
209 weights[1] * t(dβ) %*% F1 +
210 weights[2] * t(dβ2) %*% F2 +
211 weights[3] * t(dβ3) %*% F3
212 }
213 grad <- c(grad, as.double(grad_β))
214
215 grad = c(grad, #gradient on b
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 )
219
220 grad
221 },
222
223 run = function(x0)
224 {
225 "Run optimization from x0 with solver..."
226
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,
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
268 exactComp <- FALSE #TODO: global, or argument...
269
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 }