X-Git-Url: https://git.auder.net/?a=blobdiff_plain;f=pkg%2FR%2FoptimParams.R;h=c061fcf0b7d8690e901b02d8ae4943c7c940701a;hb=d08fef424150599b8095727c0f9870ca9535fb65;hp=06d16845d5784e0140fbcb4c255bef2c6196c8c6;hpb=0ad4c8de650e9f27ec3754c9cb9b2a03db5aff24;p=morpheus.git diff --git a/pkg/R/optimParams.R b/pkg/R/optimParams.R index 06d1684..c061fcf 100644 --- a/pkg/R/optimParams.R +++ b/pkg/R/optimParams.R @@ -1,17 +1,9 @@ -#' Optimize parameters -#' -#' Optimize the parameters of a mixture of logistic regressions model, possibly using -#' \code{mu <- computeMu(...)} as a partial starting point. +#' Wrapper function for OptimParams class #' #' @param K Number of populations. #' @param link The link type, 'logit' or 'probit'. -#' @param optargs a list with optional arguments: -#' \itemize{ -#' \item 'M' : list of moments of order 1,2,3: will be computed if not provided. -#' \item 'X,Y' : input/output, mandatory if moments not given -#' \item 'exact': use exact formulas when available? -#' \item weights Weights on moments when minimizing sum of squares -#' } +#' @param X Data matrix of covariables +#' @param Y Output as a binary vector #' #' @return An object 'op' of class OptimParams, initialized so that \code{op$run(x0)} #' outputs the list of optimized parameters @@ -20,8 +12,8 @@ #' \item β: regression matrix, size dxK #' \item b: intercepts, size K #' } -#' x0 is a vector containing respectively the K-1 first elements of p, then β by -#' columns, and finally b: \code{x0 = c(p[1:(K-1)],as.double(β),b)}. +#' θ0 is a vector containing respectively the K-1 first elements of p, then β by +#' columns, and finally b: \code{θ0 = c(p[1:(K-1)],as.double(β),b)}. #' #' @seealso \code{multiRun} to estimate statistics based on β, and #' \code{generateSampleIO} for I/O random generation. @@ -30,131 +22,155 @@ #' # Optimize parameters from estimated μ #' io = generateSampleIO(10000, 1/2, matrix(c(1,-2,3,1),ncol=2), c(0,0), "logit") #' μ = computeMu(io$X, io$Y, list(K=2)) -#' M <- computeMoments(io$X, io$Y) -#' o <- optimParams(2, "logit", list(M=M)) -#' x0 <- c(1/2, as.double(μ), c(0,0)) -#' par0 <- o$run(x0) +#' o <- optimParams(io$X, io$Y, 2, "logit") +#' θ0 <- list(p=1/2, β=μ, b=c(0,0)) +#' par0 <- o$run(θ0) #' # Compare with another starting point -#' x1 <- c(1/2, 2*as.double(μ), c(0,0)) -#' par1 <- o$run(x1) +#' θ1 <- list(p=1/2, β=2*μ, b=c(0,0)) +#' par1 <- o$run(θ1) #' o$f( o$linArgs(par0) ) #' o$f( o$linArgs(par1) ) #' @export -optimParams = function(K, link=c("logit","probit"), optargs=list()) +optimParams = function(X, Y, K, link=c("logit","probit")) { # Check arguments + if (!is.matrix(X) || any(is.na(X))) + stop("X: numeric matrix, no NAs") + if (!is.numeric(Y) || any(is.na(Y)) || any(Y!=0 | Y!=1)) + stop("Y: binary vector with 0 and 1 only") link <- match.arg(link) - if (!is.list(optargs)) - stop("optargs: list") - if (!is.numeric(K) || K < 2) - stop("K: integer >= 2") - - M <- optargs$M - if (is.null(M)) - { - if (is.null(optargs$X) || is.null(optargs$Y)) - stop("If moments are not provided, X and Y are required") - M <- computeMoments(optargs$X,optargs$Y) - } - - weights <- optargs$weights - if (is.null(weights)) - weights <- rep(1, 3) + if (!is.numeric(K) || K!=floor(K) || K < 2) + stop("K: integer >= 2") # Build and return optimization algorithm object - methods::new("OptimParams", "li"=link, "M1"=as.double(M[[1]]), - "M2"=as.double(M[[2]]), "M3"=as.double(M[[3]]), - "weights"=weights, "K"=as.integer(K)) + methods::new("OptimParams", "li"=link, "X"=X, + "Y"=as.integer(Y), "K"=as.integer(K)) } -# Encapsulated optimization for p (proportions), β and b (regression parameters) -# -# @field li Link, 'logit' or 'probit' -# @field M1 Estimated first-order moment -# @field M2 Estimated second-order moment (flattened) -# @field M3 Estimated third-order moment (flattened) -# @field weights Vector of moments' weights -# @field K Number of populations -# @field d Number of dimensions -# +#' Encapsulated optimization for p (proportions), β and b (regression parameters) +#' +#' Optimize the parameters of a mixture of logistic regressions model, possibly using +#' \code{mu <- computeMu(...)} as a partial starting point. +#' +#' @field li Link function, 'logit' or 'probit' +#' @field X Data matrix of covariables +#' @field Y Output as a binary vector +#' @field K Number of populations +#' @field d Number of dimensions +#' @field W Weights matrix (iteratively refined) +#' setRefClass( Class = "OptimParams", fields = list( # Inputs - li = "character", #link 'logit' or 'probit' - M1 = "numeric", #order-1 moment (vector size d) - M2 = "numeric", #M2 easier to process as a vector - M3 = "numeric", #M3 easier to process as a vector - weights = "numeric", #weights on moments + li = "character", #link function + X = "matrix", + Y = "numeric", + Mhat = "numeric", #vector of empirical moments # Dimensions K = "integer", - d = "integer" + n = "integer", + d = "integer", + # Weights matrix (generalized least square) + W = "matrix" ), methods = list( initialize = function(...) { - "Check args and initialize K, d" + "Check args and initialize K, d, W" - callSuper(...) - if (!hasArg("li") || !hasArg("M1") || !hasArg("M2") || !hasArg("M3") - || !hasArg("K")) - { + callSuper(...) + if (!hasArg("X") || !hasArg("Y") || !hasArg("K") || !hasArg("li")) stop("Missing arguments") - } + # Precompute empirical moments + M <- computeMoments(optargs$X,optargs$Y) + M1 <- as.double(M[[1]]) + M2 <- as.double(M[[2]]) + M3 <- as.double(M[[3]]) + Mhat <<- matrix(c(M1,M2,M3), ncol=1) + + n <<- nrow(X) d <<- length(M1) + W <<- diag(d+d^2+d^3) #initialize at W = Identity }, - expArgs = function(x) + expArgs = function(v) { - "Expand individual arguments from vector x" + "Expand individual arguments from vector v into a list" list( # p: dimension K-1, need to be completed - "p" = c(x[1:(K-1)], 1-sum(x[1:(K-1)])), - "β" = matrix(x[K:(K+d*K-1)], ncol=K), - "b" = x[(K+d*K):(K+(d+1)*K-1)]) + "p" = c(v[1:(K-1)], 1-sum(v[1:(K-1)])), + "β" = matrix(v[K:(K+d*K-1)], ncol=K), + "b" = v[(K+d*K):(K+(d+1)*K-1)]) }, - linArgs = function(o) + linArgs = function(L) { - " Linearize vectors+matrices into a vector x" + "Linearize vectors+matrices from list L into a vector" - c(o$p[1:(K-1)], as.double(o$β), o$b) + c(L$p[1:(K-1)], as.double(L$β), L$b) }, - f = function(x) - { - "Sum of squares (Mi - hat_Mi)^2 where Mi is obtained from formula" - - P <- expArgs(x) - p <- P$p - β <- P$β + computeW = function(θ) + { + dim <- d + d^2 + d^3 + W <<- solve( matrix( .C("Compute_Omega", + X=as.double(X), Y=as.double(Y), M=as.double(M(θ)), + pn=as.integer(n), pd=as.integer(d), + W=as.double(W), PACKAGE="morpheus")$W, nrow=dim, ncol=dim) ) + NULL #avoid returning W + }, + + M <- function(θ) + { + "Vector of moments, of size d+d^2+d^3" + + p <- θ$p + β <- θ$β λ <- sqrt(colSums(β^2)) - b <- P$b + b <- θ$b # Tensorial products β^2 = β2 and β^3 = β3 must be computed from current β1 β2 <- apply(β, 2, function(col) col %o% col) β3 <- apply(β, 2, function(col) col %o% col %o% col) - return( - weights[1] * sum( ( β %*% (p * .G(li,1,λ,b)) - M1 )^2 ) + - weights[2] * sum( ( β2 %*% (p * .G(li,2,λ,b)) - M2 )^2 ) + - weights[3] * sum( ( β3 %*% (p * .G(li,3,λ,b)) - M3 )^2 ) ) - }, + matrix(c( + β %*% (p * .G(li,1,λ,b)), + β2 %*% (p * .G(li,2,λ,b)), + β3 %*% (p * .G(li,3,λ,b))), ncol=1) + }, + + f = function(θ) + { + "Product t(Mi - hat_Mi) W (Mi - hat_Mi) with Mi(theta)" + + A <- M(θ) - Mhat + t(A) %*% W %*% A + }, - grad_f = function(x) + grad_f = function(θ) { "Gradient of f, dimension (K-1) + d*K + K = (d+2)*K - 1" - P <- expArgs(x) - p <- P$p - β <- P$β + -2 * t(grad_M(θ)) %*% getW(θ) %*% (Mhat - M(θ)) + } + + grad_M = function(θ) + { + "Gradient of the vector of moments, size (dim=)d+d^2+d^3 x K-1+K+d*K" + + L <- expArgs(θ) + p <- L$p + β <- L$β λ <- sqrt(colSums(β^2)) μ <- sweep(β, 2, λ, '/') - b <- P$b + b <- L$b + + res <- matrix(nrow=nrow(W), ncol=0) # Tensorial products β^2 = β2 and β^3 = β3 must be computed from current β1 β2 <- apply(β, 2, function(col) col %o% col) @@ -167,18 +183,14 @@ setRefClass( G4 = .G(li,4,λ,b) G5 = .G(li,5,λ,b) - # (Mi - hat_Mi)^2 ' == (Mi - hat_Mi)' 2(Mi - hat_Mi) = Mi' Fi - F1 = as.double( 2 * ( β %*% (p * G1) - M1 ) ) - F2 = as.double( 2 * ( β2 %*% (p * G2) - M2 ) ) - F3 = as.double( 2 * ( β3 %*% (p * G3) - M3 ) ) - + # Gradient on p: K-1 columns, dim rows km1 = 1:(K-1) - grad <- #gradient on p - weights[1] * t( sweep(as.matrix(β [,km1]), 2, G1[km1], '*') - G1[K] * β [,K] ) %*% F1 + - weights[2] * t( sweep(as.matrix(β2[,km1]), 2, G2[km1], '*') - G2[K] * β2[,K] ) %*% F2 + - weights[3] * t( sweep(as.matrix(β3[,km1]), 2, G3[km1], '*') - G3[K] * β3[,K] ) %*% F3 + res <- cbind(res, rbind( + t( sweep(as.matrix(β [,km1]), 2, G1[km1], '*') - G1[K] * β [,K] ), + t( sweep(as.matrix(β2[,km1]), 2, G2[km1], '*') - G2[K] * β2[,K] ), + t( sweep(as.matrix(β3[,km1]), 2, G3[km1], '*') - G3[K] * β3[,K] ))) - grad_β <- matrix(nrow=d, ncol=K) + # TODO: understand derivatives order and match the one in optim init param for (i in 1:d) { # i determines the derivated matrix dβ[2,3] @@ -205,47 +217,50 @@ setRefClass( dβ3_right[block,] <- dβ3_right[block,] + β2 dβ3 <- dβ3_left + sweep(dβ3_right, 2, p * G3, '*') - grad_β[i,] <- - weights[1] * t(dβ) %*% F1 + - weights[2] * t(dβ2) %*% F2 + - weights[3] * t(dβ3) %*% F3 + res <- cbind(res, rbind(t(dβ), t(dβ2), t(dβ3))) } - grad <- c(grad, as.double(grad_β)) - grad = c(grad, #gradient on b - weights[1] * t( sweep(β, 2, p * G2, '*') ) %*% F1 + - weights[2] * t( sweep(β2, 2, p * G3, '*') ) %*% F2 + - weights[3] * t( sweep(β3, 2, p * G4, '*') ) %*% F3 ) + # Gradient on b + res <- cbind(res, rbind( + t( sweep(β, 2, p * G2, '*') ), + t( sweep(β2, 2, p * G3, '*') ), + t( sweep(β3, 2, p * G4, '*') ))) - grad + res }, - run = function(x0) + run = function(θ0) { - "Run optimization from x0 with solver..." - - if (!is.list(x0)) - stop("x0: list") - if (is.null(x0$β)) - stop("At least x0$β must be provided") - if (!is.matrix(x0$β) || any(is.na(x0$β)) || ncol(x0$β) != K) - stop("x0$β: matrix, no NA, ncol == K") - if (is.null(x0$p)) - x0$p = rep(1/K, K-1) - else if (length(x0$p) != K-1 || sum(x0$p) > 1) - stop("x0$p should contain positive integers and sum to < 1") + "Run optimization from θ0 with solver..." + + if (!is.list(θ0)) + stop("θ0: list") + if (is.null(θ0$β)) + stop("At least θ0$β must be provided") + if (!is.matrix(θ0$β) || any(is.na(θ0$β)) || ncol(θ0$β) != K) + stop("θ0$β: matrix, no NA, ncol == K") + if (is.null(θ0$p)) + θ0$p = rep(1/K, K-1) + else if (length(θ0$p) != K-1 || sum(θ0$p) > 1) + stop("θ0$p should contain positive integers and sum to < 1") # Next test = heuristic to detect missing b (when matrix is called "beta") - if (is.null(x0$b) || all(x0$b == x0$β)) - x0$b = rep(0, K) - else if (any(is.na(x0$b))) - stop("x0$b cannot have missing values") + if (is.null(θ0$b) || all(θ0$b == θ0$β)) + θ0$b = rep(0, K) + else if (any(is.na(θ0$b))) + stop("θ0$b cannot have missing values") - op_res = constrOptim( linArgs(x0), .self$f, .self$grad_f, + op_res = constrOptim( linArgs(θ0), .self$f, .self$grad_f, ui=cbind( rbind( rep(-1,K-1), diag(K-1) ), matrix(0, nrow=K, ncol=(d+1)*K) ), ci=c(-1,rep(0,K-1)) ) + # debug: + print(computeW(expArgs(op_res$par))) + # We get a first non-trivial estimation of W + # TODO: loop, this redefine f, so that we can call constrOptim again... + # Stopping condition? N iterations? Delta <= ε ? + expArgs(op_res$par) } ) @@ -253,6 +268,7 @@ setRefClass( # Compute vectorial E[g^{(order)}(<β,x> + b)] with x~N(0,Id) (integral in R^d) # = E[g^{(order)}(z)] with z~N(b,diag(λ)) +# by numerically evaluating the integral. # # @param link Link, 'logit' or 'probit' # @param order Order of derivative @@ -264,56 +280,23 @@ setRefClass( # NOTE: weird "integral divergent" error on inputs: # link="probit"; order=2; λ=c(531.8099,586.8893,523.5816); b=c(-118.512674,-3.488020,2.109969) # Switch to pracma package for that (but it seems slow...) - - exactComp <- FALSE #TODO: global, or argument... - - if (exactComp && link == "probit") - { - # Use exact computations - sapply( seq_along(λ), function(k) { - .exactProbitIntegral(order, λ[k], b[k]) - }) - } - - else - { - # Numerical integration - sapply( seq_along(λ), function(k) { - res <- NULL - tryCatch({ - # Fast code, may fail: - res <- stats::integrate( - function(z) .deriv[[link]][[order]](λ[k]*z+b[k]) * exp(-z^2/2) / sqrt(2*pi), - lower=-Inf, upper=Inf )$value - }, error = function(e) { - # Robust slow code, no fails observed: - sink("/dev/null") #pracma package has some useless printed outputs... - res <- pracma::integral( - function(z) .deriv[[link]][[order]](λ[k]*z+b[k]) * exp(-z^2/2) / sqrt(2*pi), - xmin=-Inf, xmax=Inf, method="Kronrod") - sink() - }) - res - }) - } -} - -# TODO: check these computations (wrong atm) -.exactProbitIntegral <- function(order, λ, b) -{ - c1 = (1/sqrt(2*pi)) * exp( -.5 * b/((λ^2+1)^2) ) - if (order == 1) - return (c1) - c2 = b - λ^2 / (λ^2+1) - if (order == 2) - return (c1 * c2) - if (order == 3) - return (c1 * (λ^2 - 1 + c2^2)) - if (order == 4) - return ( (c1*c2/((λ^2+1)^2)) * (-λ^4*((b+1)^2+1) - - 2*λ^3 + λ^2*(2-2*b*(b-1)) + 6*λ + 3 - b^2) ) - if (order == 5) #only remaining case... - return ( c1 * (3*λ^4+c2^4+6*c1^2*(λ^2-1) - 6*λ^2 + 6) ) + sapply( seq_along(λ), function(k) { + res <- NULL + tryCatch({ + # Fast code, may fail: + res <- stats::integrate( + function(z) .deriv[[link]][[order]](λ[k]*z+b[k]) * exp(-z^2/2) / sqrt(2*pi), + lower=-Inf, upper=Inf )$value + }, error = function(e) { + # Robust slow code, no fails observed: + sink("/dev/null") #pracma package has some useless printed outputs... + res <- pracma::integral( + function(z) .deriv[[link]][[order]](λ[k]*z+b[k]) * exp(-z^2/2) / sqrt(2*pi), + xmin=-Inf, xmax=Inf, method="Kronrod") + sink() + }) + res + }) } # Derivatives list: g^(k)(x) for links 'logit' and 'probit' @@ -321,7 +304,7 @@ setRefClass( .deriv <- list( "probit"=list( # 'probit' derivatives list; - # TODO: exact values for the integral E[g^(k)(λz+b)] + # NOTE: exact values for the integral E[g^(k)(λz+b)] could be computed function(x) exp(-x^2/2)/(sqrt(2*pi)), #g' function(x) exp(-x^2/2)/(sqrt(2*pi)) * -x, #g'' function(x) exp(-x^2/2)/(sqrt(2*pi)) * ( x^2 - 1), #g^(3)