Remove weights for now; planning generalization with matrix W
[morpheus.git] / pkg / R / optimParams.R
... / ...
CommitLineData
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
43optimParams = 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 # Build and return optimization algorithm object
61 methods::new("OptimParams", "li"=link, "M1"=as.double(M[[1]]),
62 "M2"=as.double(M[[2]]), "M3"=as.double(M[[3]]), "K"=as.integer(K))
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#
74setRefClass(
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",
85 d = "integer",
86 # Weights matrix (generalized least square)
87 W = "matrix"
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)
103 W <<- diag(d+d^2+d^3) #initialize at W = Identity
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 {
126 "Product t(Mi - hat_Mi) W (Mi - hat_Mi) with Mi(theta)"
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(
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 ) )
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
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
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
204 grad_β[i,] <- t(dβ) %*% F1 + t(dβ2) %*% F2 + t(dβ3) %*% F3
205 }
206 grad <- c(grad, as.double(grad_β))
207
208 grad = c(grad, #gradient on b
209 t( sweep(β, 2, p * G2, '*') ) %*% F1 +
210 t( sweep(β2, 2, p * G3, '*') ) %*% F2 +
211 t( sweep(β3, 2, p * G4, '*') ) %*% F3 )
212
213 grad
214 },
215
216 run = function(x0)
217 {
218 "Run optimization from x0 with solver..."
219
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,
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
261 exactComp <- FALSE #TODO: global, or argument...
262
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}