--- /dev/null
+*.pdf filter=fat
+*.html filter=fat
+*.zip filter=fat
+*.tar.xz filter=fat
+*.data filter=fat
+*.png filter=fat
--- /dev/null
+[rsync]
+remote = gitfat@auder.net:~/files/morpheus
--- /dev/null
+#ignore roxygen2 generated files
+NAMESPACE
+*.Rd
+!*-package.Rd
+
+#ignore temporary files
+*~
+*.swp
+
+#ignore package for CRAN
+/pkg-cran/
+
+#ignore R session files
+.Rhistory
+.RData
+
+#ignore R CMD build/check genrated files
+/*.Rcheck/
+/*.tar.gz
+
+#ignore generated object files
+*.[oa]
+*.so
--- /dev/null
+# Estimate parameters of mixtures of logistic regressions
+
+This code is one applied part in the PhD thesis of [Mor-Absa Loum](https://fr.linkedin.com/in/mor-absa-loum-3372aa73)
+
+## Description
+
+Mixture of lOgistic Regressions Parameters (H)Estimation with (U)Spectral methods.
+The main methods take d-dimensional inputs + a vector of binary outputs, and return
+parameters according to the GLMs mixture model (please see the package vignette).
+
+NOTE: greek unicode letters are used in the code, because it's much nicer to write λ, β
+and Σ than lambda, beta and Sigma - and also more importantly because it works well :)
+...However CRAN demands ASCII files - therefore the presence of to-cran.sh script.
+
+## Example
+
+Install the package, then ?multiRun is a possible starting point.
--- /dev/null
+simus point initial chaque M-C --> 3 ou 5 ? multiples de mu ? perturbés ?
+n petit, mu = mauvais point initial ? plusieurs ? combien ? recommandations ?
+n grand : mu OK ? pour quel n en fonction de d ?
+
+Pour la vérification de M_3 indiquée par Francis, juste un point:
+Si on note B la matrice qui contient les beta en colonne, et si on note A l'inverse de B, c'est
+A M_3 (\rho_l) A^{T} qui est diagonale
+(car M_3 (\rho_l) = B D_l B^{T} pour une matrice D_l diagonale).
+
+Partir du vrai beta + vrais param :: voir si on converge loin de beta (n diminue...)
+--> partir du vrai mu, aussi.
+betas sparses, bcp de 0 ?! (comment ça se comporte ?)
--- /dev/null
+#$# git-fat 16f72d1107b66f90f27dd19ae241fc1a528971a5 478827
--- /dev/null
+FLXMRglm <- function(formula=.~., family=gaussian, offset=NULL)
+{
+ if (is.character(family))
+ family <- get(family, mode = "function", envir = parent.frame())
+ if (is.function(family))
+ family <- family()
+ if (is.null(family$family)) {
+ print(family)
+ stop("'family' not recognized")
+ }
+ glmrefit <- function(x, y, w) {
+ fit <- c(glm.fit(x, y, weights=w, offset=offset,
+ family=family),
+ list(call = sys.call(), offset = offset,
+ control = eval(formals(glm.fit)$control),
+ method = "weighted.glm.fit"))
+ fit$df.null <- sum(w) + fit$df.null - fit$df.residual - fit$rank
+ fit$df.residual <- sum(w) - fit$rank
+ fit$x <- x
+ fit
+ }
+
+ z <- new("FLXMRglm", weighted=TRUE, formula=formula,
+ name=paste("FLXMRglm", family$family, sep=":"), offset = offset,
+ family=family$family, refit=glmrefit)
+ z@preproc.y <- function(x){
+ if (ncol(x) > 1)
+ stop(paste("for the", family$family, "family y must be univariate"))
+ x
+ }
+
+ if(family$family=="gaussian"){
+ z@defineComponent <- function(para) {
+ predict <- function(x, ...) {
+ dotarg = list(...)
+ if("offset" %in% names(dotarg)) offset <- dotarg$offset
+ p <- x %*% para$coef
+ if (!is.null(offset)) p <- p + offset
+ family$linkinv(p)
+ }
+
+ logLik <- function(x, y, ...)
+ dnorm(y, mean=predict(x, ...), sd=para$sigma, log=TRUE)
+
+ new("FLXcomponent",
+ parameters=list(coef=para$coef, sigma=para$sigma),
+ logLik=logLik, predict=predict,
+ df=para$df)
+ }
+
+ z@fit <- function(x, y, w, component){
+ fit <- glm.fit(x, y, w=w, offset=offset, family = family)
+ z@defineComponent(para = list(coef = coef(fit), df = ncol(x)+1,
+ sigma = sqrt(sum(fit$weights * fit$residuals^2 /
+ mean(fit$weights))/ (nrow(x)-fit$rank))))
+ }
+ }
+ else if(family$family=="binomial"){
+ z@preproc.y <- function(x){
+ if (ncol(x) != 2)
+ stop("for the binomial family, y must be a 2 column matrix\n",
+ "where col 1 is no. successes and col 2 is no. failures")
+ if (any(x < 0))
+ stop("negative values are not allowed for the binomial family")
+ x
+ }
+ z@defineComponent <- function(para) {
+ predict <- function(x, ...) {
+ dotarg = list(...)
+ if("offset" %in% names(dotarg)) offset <- dotarg$offset
+ p <- x %*% para$coef
+ if (!is.null(offset)) p <- p + offset
+ family$linkinv(p)
+ }
+ logLik <- function(x, y, ...)
+ dbinom(y[,1], size=rowSums(y), prob=predict(x, ...), log=TRUE)
+
+ new("FLXcomponent",
+ parameters=list(coef=para$coef),
+ logLik=logLik, predict=predict,
+ df=para$df)
+ }
+
+ z@fit <- function(x, y, w, component){
+ fit <- glm.fit(x, y, weights=w, family=family, offset=offset, start=component$coef)
+ z@defineComponent(para = list(coef = coef(fit), df = ncol(x)))
+ }
+ }
+ else if(family$family=="poisson"){
+ z@defineComponent <- function(para) {
+ predict <- function(x, ...) {
+ dotarg = list(...)
+ if("offset" %in% names(dotarg)) offset <- dotarg$offset
+ p <- x %*% para$coef
+ if (!is.null(offset)) p <- p + offset
+ family$linkinv(p)
+ }
+ logLik <- function(x, y, ...)
+ dpois(y, lambda=predict(x, ...), log=TRUE)
+
+ new("FLXcomponent",
+ parameters=list(coef=para$coef),
+ logLik=logLik, predict=predict,
+ df=para$df)
+ }
+
+ z@fit <- function(x, y, w, component){
+ fit <- glm.fit(x, y, weights=w, family=family, offset=offset, start=component$coef)
+ z@defineComponent(para = list(coef = coef(fit), df = ncol(x)))
+ }
+ }
+ else if(family$family=="Gamma"){
+ z@defineComponent <- function(para) {
+ predict <- function(x, ...) {
+ dotarg = list(...)
+ if("offset" %in% names(dotarg)) offset <- dotarg$offset
+ p <- x %*% para$coef
+ if (!is.null(offset)) p <- p + offset
+ family$linkinv(p)
+ }
+ logLik <- function(x, y, ...)
+ dgamma(y, shape = para$shape, scale=predict(x, ...)/para$shape, log=TRUE)
+
+ new("FLXcomponent",
+ parameters = list(coef = para$coef, shape = para$shape),
+ predict = predict, logLik = logLik,
+ df = para$df)
+ }
+
+ z@fit <- function(x, y, w, component){
+ fit <- glm.fit(x, y, weights=w, family=family, offset=offset, start=component$coef)
+ z@defineComponent(para = list(coef = coef(fit), df = ncol(x)+1,
+ shape = sum(fit$prior.weights)/fit$deviance))
+
+ }
+ }
+ else stop(paste("Unknown family", family))
+ z
+}
--- /dev/null
+library("flexmix")
+data("tribolium", package = "flexmix")
+set.seed(1234)
+## uses FLXMRglm() from package.
+tribMix <- initFlexmix(cbind(Remaining, Total - Remaining) ~ Species,
+ k = 2, nrep = 5, data = tribolium,
+ model = FLXMRglm(family = "binomial"))
+parameters(tribMix); logLik(tribMix)
+source("FLXMRglm.R")
+set.seed(1234)
+## uses FLXMRglm() from source file which allows to specify link.
+tribMixNew <- initFlexmix(cbind(Remaining, Total - Remaining) ~ Species,
+ k = 2, nrep = 5, data = tribolium,
+ model = FLXMRglm(family = "binomial"))
+parameters(tribMixNew); logLik(tribMixNew)
+set.seed(1234)
+tribMixlogit <- initFlexmix(cbind(Remaining, Total - Remaining) ~ Species,
+ k = 2, nrep = 5, data = tribolium,
+ model = FLXMRglm(family = binomial(link = logit)))
+parameters(tribMixlogit); logLik(tribMixlogit)
+set.seed(1234)
+tribMixprobit <- initFlexmix(cbind(Remaining, Total - Remaining) ~ Species,
+ k = 2, nrep = 5, data = tribolium,
+ model = FLXMRglm(family = binomial(link = probit)))
+parameters(tribMixprobit); logLik(tribMixprobit)
--- /dev/null
+Package: morpheus
+Title: Estimate Parameters of Mixtures of Logistic Regressions
+Description: Mixture of lOgistic Regressions Parameters (H)Estimation with
+ (U)Spectral methods. The main methods take d-dimensional inputs and a vector
+ of binary outputs, and return parameters according to the GLMs mixture model
+ (please refer to the package vignette).
+Version: 0.2-0
+Author: Benjamin Auder <Benjamin.Auder@math.u-psud.fr> [aut,cre],
+ Mor-Absa Loum <Mor-Absa.Loum@math.u-psud.fr> [aut]
+Maintainer: Benjamin Auder <Benjamin.Auder@math.u-psud.fr>
+Depends:
+ R (>= 3.0.0),
+Imports:
+ MASS,
+ jointDiag,
+ methods,
+ pracma
+Suggests:
+ flexmix,
+ parallel,
+ testthat,
+ roxygen2,
+ tensor,
+ nloptr
+License: MIT + file LICENSE
+RoxygenNote: 5.0.1
+Collate:
+ 'utils.R'
+ 'A_NAMESPACE.R'
+ 'computeMu.R'
+ 'multiRun.R'
+ 'optimParams.R'
+ 'plot.R'
+ 'sampleIO.R'
--- /dev/null
+Copyright (c) 2016-2017, Mor-Absa Loum ; 2017, Benjamin Auder
+
+Permission is hereby granted, free of charge, to any person obtaining
+a copy of this software and associated documentation files (the
+"Software"), to deal in the Software without restriction, including
+without limitation the rights to use, copy, modify, merge, publish,
+distribute, sublicense, and/or sell copies of the Software, and to
+permit persons to whom the Software is furnished to do so, subject to
+the following conditions:
+
+The above copyright notice and this permission notice shall be
+included in all copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
+MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
+NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE
+LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION
+OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
+WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
--- /dev/null
+#' @include utils.R
+#'
+#' @useDynLib morpheus
+#'
+#' @importFrom jointDiag ajd
+#' @importFrom stats rbinom rmultinom rnorm pnorm runif integrate
+#' @importFrom graphics boxplot barplot hist par
+#' @importFrom methods new
+#' @importFrom pracma integral
+#'
+NULL
--- /dev/null
+#' Compute μ
+#'
+#' Estimate the normalized columns μ of the β matrix parameter in a mixture of
+#' logistic regressions models, with a spectral method described in the package vignette.
+#'
+#' @param X Matrix of input data (size nxd)
+#' @param Y Vector of binary outputs (size n)
+#' @param optargs List of optional argument:
+#' \itemize{
+#' \item 'jd_method', joint diagonalization method from the package jointDiag:
+#' 'uwedge' (default) or 'jedi'.
+#' \item 'jd_nvects', number of random vectors for joint-diagonalization
+#' (or 0 for p=d, canonical basis by default)
+#' \item 'M', moments of order 1,2,3: will be computed if not provided.
+#' \item 'K', number of populations (estimated with ranks of M2 if not given)
+#' }
+#'
+#' @return The estimated normalized parameters as columns of a matrix μ of size dxK
+#'
+#' @seealso \code{multiRun} to estimate statistics based on μ,
+#' and \code{generateSampleIO} for I/O random generation.
+#'
+#' @examples
+#' io = generateSampleIO(10000, 1/2, matrix(c(1,0,0,1),ncol=2), c(0,0), "probit")
+#' μ = computeMu(io$X, io$Y, list(K=2)) #or just X and Y for estimated K
+#' @export
+computeMu = function(X, Y, optargs=list())
+{
+ if (!is.matrix(X) || !is.numeric(X) || any(is.na(X)))
+ stop("X: real matrix, no NA")
+ n = nrow(X)
+ d = ncol(X)
+ if (!is.numeric(Y) || length(Y)!=n || any(Y!=0 & Y!=1))
+ stop("Y: vector of 0 and 1, size nrow(X), no NA")
+ if (!is.list(optargs))
+ stop("optargs: list")
+
+ # Step 0: Obtain the empirically estimated moments tensor, estimate also K
+ M = if (is.null(optargs$M)) computeMoments(X,Y) else optargs$M
+ K = optargs$K
+ if (is.null(K))
+ {
+ # TODO: improve this basic heuristic
+ Σ = svd(M[[2]])$d
+ large_ratio <- ( abs(Σ[-d] / Σ[-1]) > 3 )
+ K <- if (any(large_ratio)) max(2, which.min(large_ratio)) else d
+ }
+
+ # Step 1: generate a family of d matrices to joint-diagonalize to increase robustness
+ d = ncol(X)
+ fixed_design = FALSE
+ jd_nvects = ifelse(!is.null(optargs$jd_nvects), optargs$jd_nvects, 0)
+ if (jd_nvects == 0)
+ {
+ jd_nvects = d
+ fixed_design = TRUE
+ }
+ M2_t = array(dim=c(d,d,jd_nvects))
+ for (i in seq_len(jd_nvects))
+ {
+ rho = if (fixed_design) c(rep(0,i-1),1,rep(0,d-i)) else normalize( rnorm(d) )
+ M2_t[,,i] = .T_I_I_w(M[[3]],rho)
+ }
+
+ # Step 2: obtain factors u_i (and their inverse) from the joint diagonalisation of M2_t
+ jd_method = ifelse(!is.null(optargs$jd_method), optargs$jd_method, "uwedge")
+ V =
+ if (jd_nvects > 1) {
+ #NOTE: increasing itermax does not help to converge, thus we suppress warnings
+ suppressWarnings({jd = jointDiag::ajd(M2_t, method=jd_method)})
+# if (jd_method=="uwedge") jd$B else solve(jd$A)
+ if (jd_method=="uwedge") jd$B else MASS::ginv(jd$A)
+ }
+ else
+ eigen(M2_t[,,1])$vectors
+
+ # Step 3: obtain final factors from joint diagonalisation of T(I,I,u_i)
+ M2_t = array(dim=c(d,d,K))
+ for (i in seq_len(K))
+ M2_t[,,i] = .T_I_I_w(M[[3]],V[,i])
+ suppressWarnings({jd = jointDiag::ajd(M2_t, method=jd_method)})
+# U = if (jd_method=="uwedge") solve(jd$B) else jd$A
+ U = if (jd_method=="uwedge") MASS::ginv(jd$B) else jd$A
+ μ = normalize(U[,1:K])
+
+ # M1 also writes M1 = sum_k coeff_k * μ_k, where coeff_k >= 0
+ # ==> search decomposition of vector M1 onto the (truncated) basis μ (of size dxK)
+ # This is a linear system μ %*% C = M1 with C of size K ==> C = psinv(μ) %*% M1
+ C = MASS::ginv(μ) %*% M[[1]]
+ μ[,C < 0] = - μ[,C < 0]
+ μ
+}
--- /dev/null
+#' multiRun
+#'
+#' Estimate N times some parameters, outputs of some list of functions. This method is
+#' thus very generic, allowing typically bootstrap or Monte-Carlo estimations of matrices
+#' μ or β. Passing a list of functions opens the possibility to compare them on a fair
+#' basis (exact same inputs). It's even possible to compare methods on some deterministic
+#' design of experiments.
+#'
+#' @param fargs List of arguments for the estimation functions
+#' @param estimParams List of nf function(s) to apply on fargs - shared signature
+#' @param prepareArgs Prepare arguments for the functions inside estimParams
+#' @param N Number of runs
+#' @param ncores Number of cores for parallel runs (<=1: sequential)
+#' @param verbose TRUE to indicate runs + methods numbers
+#'
+#' @return A list of nf aggregates of N results (matrices).
+#'
+#' @examples
+#' \dontrun{
+#' β <- matrix(c(1,-2,3,1),ncol=2)
+#'
+#' # Bootstrap + computeMu, morpheus VS flexmix ; assumes fargs first 3 elts X,Y,K
+#' io <- generateSampleIO(n=1000, p=1/2, β=β, b=c(0,0), "logit")
+#' μ <- normalize(β)
+#' res <- multiRun(list(X=io$X,Y=io$Y,optargs=list(K=2,jd_nvects=0)), list(
+#' # morpheus
+#' function(fargs) {
+#' library(morpheus)
+#' ind <- fargs$ind
+#' computeMu(fargs$X[ind,],fargs$Y[ind],fargs$optargs)
+#' },
+#' # flexmix
+#' function(fargs) {
+#' library(flexmix)
+#' ind <- fargs$ind
+#' K <- fargs$optargs$K
+#' dat = as.data.frame( cbind(fargs$Y[ind],fargs$X[ind,]) )
+#' out = refit( flexmix( cbind(V1, 1 - V1) ~ 0+., data=dat, k=K,
+#' model=FLXMRglm(family="binomial") ) )
+#' normalize( matrix(out@@coef[1:(ncol(fargs$X)*K)], ncol=K) )
+#' } ),
+#' prepareArgs = function(fargs) {
+#' fargs$ind <- sample(1:nrow(fargs$X),replace=TRUE)
+#' fargs
+#' }, N=10, ncores=3)
+#' for (i in 1:2)
+#' res[[i]] <- alignMatrices(res[[i]], ref=μ, ls_mode="exact")
+#'
+#' # Monte-Carlo + optimParams from X,Y, morpheus VS flexmix ; first args n,p,β,b
+#' res <- multiRun(list(n=1000,p=1/2,β=β,b=c(0,0),optargs=list(link="logit")),list(
+#' # morpheus
+#' function(fargs) {
+#' library(morpheus)
+#' K <- fargs$optargs$K
+#' μ <- computeMu(fargs$X, fargs$Y, fargs$optargs)
+#' V <- list( p=rep(1/K,K-1), β=μ, b=c(0,0) )
+#' optimParams(V,fargs$optargs)$β
+#' },
+#' # flexmix
+#' function(fargs) {
+#' library(flexmix)
+#' K <- fargs$optargs$K
+#' dat <- as.data.frame( cbind(fargs$Y,fargs$X) )
+#' out <- refit( flexmix( cbind(V1, 1 - V1) ~ 0+., data=dat, k=K,
+#' model=FLXMRglm(family="binomial") ) )
+#' sapply( seq_len(K), function(i) as.double( out@@components[[1]][[i]][,1] ) )
+#' } ),
+#' prepareArgs = function(fargs) {
+#' library(morpheus)
+#' io = generateSampleIO(fargs$n, fargs$p, fargs$β, fargs$b, fargs$optargs$link)
+#' fargs$X = io$X
+#' fargs$Y = io$Y
+#' fargs$optargs$K = ncol(fargs$β)
+#' fargs$optargs$M = computeMoments(io$X,io$Y)
+#' fargs
+#' }, N=10, ncores=3)
+#' for (i in 1:2)
+#' res[[i]] <- alignMatrices(res[[i]], ref=β, ls_mode="exact")}
+#' @export
+multiRun <- function(fargs, estimParams,
+ prepareArgs = function(x) x, N=10, ncores=3, agg=lapply, verbose=FALSE)
+{
+ if (!is.list(fargs))
+ stop("fargs: list")
+ # No checks on fargs: supposedly done in estimParams[[i]]()
+ if (!is.list(estimParams))
+ estimParams = list(estimParams)
+ # Verify that the provided parameters estimations are indeed functions
+ lapply(seq_along(estimParams), function(i) {
+ if (!is.function(estimParams[[i]]))
+ stop("estimParams: list of function(fargs)")
+ })
+ if (!is.numeric(N) || N < 1)
+ stop("N: positive integer")
+
+ estimParamAtIndex <- function(index)
+ {
+ fargs <- prepareArgs(fargs)
+ if (verbose)
+ cat("Run ",index,"\n")
+ lapply(seq_along(estimParams), function(i) {
+ if (verbose)
+ cat(" Method ",i,"\n")
+ out <- estimParams[[i]](fargs)
+ if (is.list(out))
+ do.call(rbind, out)
+ else
+ out
+ })
+ }
+
+ if (ncores > 1)
+ {
+ cl = parallel::makeCluster(ncores, outfile="")
+ parallel::clusterExport(cl, c("fargs","verbose"), environment())
+ list_res = parallel::clusterApplyLB(cl, 1:N, estimParamAtIndex)
+ parallel::stopCluster(cl)
+ }
+ else
+ list_res = lapply(1:N, estimParamAtIndex)
+
+ # De-interlace results: output one list per function
+ nf <- length(estimParams)
+ lapply( seq_len(nf), function(i) lapply(seq_len(N), function(j) list_res[[j]][[i]]) )
+}
--- /dev/null
+#' Optimize parameters
+#'
+#' Optimize the parameters of a mixture of logistic regressions model, possibly using
+#' \code{mu <- computeMu(...)} as a partial starting point.
+#'
+#' @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?
+#' }
+#'
+#' @return An object 'op' of class OptimParams, initialized so that \code{op$run(x0)}
+#' outputs the list of optimized parameters
+#' \itemize{
+#' \item p: proportions, size K
+#' \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)}.
+#'
+#' @seealso \code{multiRun} to estimate statistics based on β, and
+#' \code{generateSampleIO} for I/O random generation.
+#'
+#' @examples
+#' # 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)
+#' # Compare with another starting point
+#' x1 <- c(1/2, 2*as.double(μ), c(0,0))
+#' par1 <- o$run(x1)
+#' o$f( o$linArgs(par0) )
+#' o$f( o$linArgs(par1) )
+#' @export
+optimParams = function(K, link=c("logit","probit"), optargs=list())
+{
+ # Check arguments
+ 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)
+ }
+
+ # TODO: field?!
+ exactComp <<- optargs$exact
+ if (is.null(exactComp))
+ exactComp <<- FALSE
+
+ # 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]]), "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 K Number of populations
+# @field d Number of dimensions
+#
+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
+ # Dimensions
+ K = "integer",
+ d = "integer"
+ ),
+
+ methods = list(
+ initialize = function(...)
+ {
+ "Check args and initialize K, d"
+
+ callSuper(...)
+ if (!hasArg("li") || !hasArg("M1") || !hasArg("M2") || !hasArg("M3")
+ || !hasArg("K"))
+ {
+ stop("Missing arguments")
+ }
+
+ d <<- length(M1)
+ },
+
+ expArgs = function(x)
+ {
+ "Expand individual arguments from vector x"
+
+ 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)])
+ },
+
+ linArgs = function(o)
+ {
+ " Linearize vectors+matrices into a vector x"
+
+ c(o$p[1:(K-1)], as.double(o$β), o$b)
+ },
+
+ f = function(x)
+ {
+ "Sum of squares (Mi - hat_Mi)^2 where Mi is obtained from formula"
+
+ P <- expArgs(x)
+ p <- P$p
+ β <- P$β
+ λ <- sqrt(colSums(β^2))
+ b <- P$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(
+ sum( ( β %*% (p * .G(li,1,λ,b)) - M1 )^2 ) +
+ sum( ( β2 %*% (p * .G(li,2,λ,b)) - M2 )^2 ) +
+ sum( ( β3 %*% (p * .G(li,3,λ,b)) - M3 )^2 ) )
+ },
+
+ grad_f = function(x)
+ {
+ "Gradient of f, dimension (K-1) + d*K + K = (d+2)*K - 1"
+
+ P <- expArgs(x)
+ p <- P$p
+ β <- P$β
+ λ <- sqrt(colSums(β^2))
+ μ <- sweep(β, 2, λ, '/')
+ b <- P$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)
+
+ # Some precomputations
+ G1 = .G(li,1,λ,b)
+ G2 = .G(li,2,λ,b)
+ G3 = .G(li,3,λ,b)
+ 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 ) )
+
+ km1 = 1:(K-1)
+ grad <- #gradient on p
+ t( sweep(as.matrix(β [,km1]), 2, G1[km1], '*') - G1[K] * β [,K] ) %*% F1 +
+ t( sweep(as.matrix(β2[,km1]), 2, G2[km1], '*') - G2[K] * β2[,K] ) %*% F2 +
+ t( sweep(as.matrix(β3[,km1]), 2, G3[km1], '*') - G3[K] * β3[,K] ) %*% F3
+
+ grad_β <- matrix(nrow=d, ncol=K)
+ for (i in 1:d)
+ {
+ # i determines the derivated matrix dβ[2,3]
+
+ dβ_left <- sweep(β, 2, p * G3 * β[i,], '*')
+ dβ_right <- matrix(0, nrow=d, ncol=K)
+ block <- i
+ dβ_right[block,] <- dβ_right[block,] + 1
+ dβ <- dβ_left + sweep(dβ_right, 2, p * G1, '*')
+
+ dβ2_left <- sweep(β2, 2, p * G4 * β[i,], '*')
+ dβ2_right <- do.call( rbind, lapply(1:d, function(j) {
+ sweep(dβ_right, 2, β[j,], '*')
+ }) )
+ block <- ((i-1)*d+1):(i*d)
+ dβ2_right[block,] <- dβ2_right[block,] + β
+ dβ2 <- dβ2_left + sweep(dβ2_right, 2, p * G2, '*')
+
+ dβ3_left <- sweep(β3, 2, p * G5 * β[i,], '*')
+ dβ3_right <- do.call( rbind, lapply(1:d, function(j) {
+ sweep(dβ2_right, 2, β[j,], '*')
+ }) )
+ block <- ((i-1)*d*d+1):(i*d*d)
+ dβ3_right[block,] <- dβ3_right[block,] + β2
+ dβ3 <- dβ3_left + sweep(dβ3_right, 2, p * G3, '*')
+
+ grad_β[i,] <- t(dβ) %*% F1 + t(dβ2) %*% F2 + t(dβ3) %*% F3
+ }
+ grad <- c(grad, as.double(grad_β))
+
+ grad = c(grad, #gradient on b
+ t( sweep(β, 2, p * G2, '*') ) %*% F1 +
+ t( sweep(β2, 2, p * G3, '*') ) %*% F2 +
+ t( sweep(β3, 2, p * G4, '*') ) %*% F3 )
+
+ grad
+ },
+
+ run = function(x0)
+ {
+ "Run optimization from x0 with solver..."
+
+ if (!is.numeric(x0) || any(is.na(x0)) || length(x0) != (d+2)*K-1
+ || any(x0[1:(K-1)] < 0) || sum(x0[1:(K-1)]) > 1)
+ {
+ stop("x0: numeric vector, no NA, length (d+2)*K-1, sum(x0[1:(K-1) >= 0]) <= 1")
+ }
+
+ op_res = constrOptim( x0, .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)) )
+
+ expArgs(op_res$par)
+ }
+ )
+)
+
+# 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(λ))
+#
+# @param link Link, 'logit' or 'probit'
+# @param order Order of derivative
+# @param λ Norm of columns of β
+# @param b Intercept
+#
+.G <- function(link, order, λ, b)
+{
+ # 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...)
+
+ 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) )
+}
+
+# Derivatives list: g^(k)(x) for links 'logit' and 'probit'
+#
+.deriv <- list(
+ "probit"=list(
+ # 'probit' derivatives list;
+ # TODO: exact values for the integral E[g^(k)(λz+b)]
+ 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)
+ function(x) exp(-x^2/2)/(sqrt(2*pi)) * (-x^3 + 3*x), #g^(4)
+ function(x) exp(-x^2/2)/(sqrt(2*pi)) * ( x^4 - 6*x^2 + 3) #g^(5)
+ ),
+ "logit"=list(
+ # Sigmoid derivatives list, obtained with http://www.derivative-calculator.net/
+ # @seealso http://www.ece.uc.edu/~aminai/papers/minai_sigmoids_NN93.pdf
+ function(x) {e=exp(x); .zin(e /(e+1)^2)}, #g'
+ function(x) {e=exp(x); .zin(e*(-e + 1) /(e+1)^3)}, #g''
+ function(x) {e=exp(x); .zin(e*( e^2 - 4*e + 1) /(e+1)^4)}, #g^(3)
+ function(x) {e=exp(x); .zin(e*(-e^3 + 11*e^2 - 11*e + 1) /(e+1)^5)}, #g^(4)
+ function(x) {e=exp(x); .zin(e*( e^4 - 26*e^3 + 66*e^2 - 26*e + 1)/(e+1)^6)} #g^(5)
+ )
+)
+
+# Utility for integration: "[return] zero if [argument is] NaN" (Inf / Inf divs)
+#
+# @param x Ratio of polynoms of exponentials, as in .S[[i]]
+#
+.zin <- function(x)
+{
+ x[is.nan(x)] <- 0.
+ x
+}
--- /dev/null
+# extractParam
+#
+# Extract successive values of a projection of the parameter(s)
+#
+# @inheritParams plotHist
+#
+extractParam <- function(mr, x=1, y=1)
+{
+ # Obtain L vectors where L = number of res lists in mr
+ lapply( mr, function(mr_list) {
+ sapply(mr_list, function(m) m[x,y])
+ } )
+}
+
+#' plotHist
+#'
+#' Plot histogram
+#'
+#' @param mr Output of multiRun(), list of lists of functions results
+#' @param x Row index of the element inside the aggregated parameter
+#' @param y Colomn index of the element inside the aggregated parameter
+#'
+#' @examples
+#' \dontrun{
+#' β <- matrix(c(1,-2,3,1),ncol=2)
+#' mr <- multiRun(...) #see bootstrap example in ?multiRun : return lists of mu_hat
+#' μ <- normalize(β)
+#' for (i in 1:2)
+#' mr[[i]] <- alignMatrices(res[[i]], ref=μ, ls_mode="exact")
+#' plotHist(mr, 2, 1) #second row, first column}
+#' @export
+plotHist <- function(mr, x, y)
+{
+ params <- extractParam(mr, x, y)
+ L = length(params)
+ # Plot histograms side by side
+ par(mfrow=c(1,L), cex.axis=1.5, cex.lab=1.5, mar=c(4.7,5,1,1))
+ for (i in 1:L)
+ hist(params[[i]], breaks=40, freq=FALSE, xlab="Parameter value", ylab="Density")
+}
+
+#' plotBox
+#'
+#' Draw boxplot
+#'
+#' @inheritParams plotHist
+#'
+#' @examples
+#' #See example in ?plotHist
+#' @export
+plotBox <- function(mr, x, y)
+{
+ params <- extractParam(mr, x, y)
+ L = length(params)
+ # Plot boxplots side by side
+ par(mfrow=c(1,L), cex.axis=1.5, cex.lab=1.5, mar=c(4.7,5,1,1))
+ for (i in 1:L)
+ boxplot(params[[i]], ylab="Parameter value")
+}
+
+#' plotCoefs
+#'
+#' Draw coefs estimations + standard deviations
+#'
+#' @inheritParams plotHist
+#' @param params True value of parameters matrix
+#'
+#' @examples
+#' #See example in ?plotHist
+#' @export
+plotCoefs <- function(mr, params)
+{
+ nf <- length(mr)
+ L <- nrow(mr[[1]][[1]])
+ K <- ncol(mr[[1]][[1]])
+
+ params_hat <- vector("list", nf)
+ stdev <- vector("list", nf)
+ for (i in 1:nf)
+ {
+ params_hat[[i]] <- matrix(nrow=L, ncol=K)
+ stdev[[i]] <- matrix(nrow=L, ncol=K)
+ }
+ for (x in 1:L)
+ {
+ for (y in 1:K)
+ {
+ estims <- extractParam(mr, x, y)
+ for (i in 1:nf)
+ {
+ params_hat[[i]][x,y] <- mean(estims[[i]])
+# stdev[[i]][x,y] <- sqrt( mean( (estims[[i]] - params[x,y])^2 ) )
+ # HACK remove extreme quantile in estims[[i]] before computing sd()
+ stdev[[i]][x,y] <- sd( estims[[i]] [ estims[[i]] < max(estims[[i]]) & estims[[i]] > min(estims[[i]]) ] )
+ }
+ }
+ }
+
+ par(mfrow=c(1,nf), cex.axis=1.5, cex.lab=1.5, mar=c(4.7,5,1,1))
+ params <- as.double(params)
+ o <- order(params)
+ for (i in 1:nf)
+ {
+ avg_param <- as.double(params_hat[[i]])
+ std_param <- as.double(stdev[[i]])
+ matplot(cbind(params[o],avg_param[o],avg_param[o]+std_param[o],avg_param[o]-std_param[o]),
+ col=c(2,1,1,1), lty=c(1,1,2,2), type="l", lwd=2, xlab="param", ylab="value")
+ }
+
+ #print(o) #not returning o to avoid weird Jupyter issue... (TODO:)
+}
+
+#' plotQn
+#'
+#' Draw 3D map of objective function values
+#'
+#' @param N Number of starting points
+#' @param β Regression matrix (target)
+#' @param link Link function (logit or probit)
+#'
+#' @export
+plotQn <- function(N, n, p, β, b, link)
+{
+ d <- nrow(β)
+ K <- ncol(β)
+ io <- generateSampleIO(n, p, β, b, link)
+ op <- optimParams(K, link, list(X=io$X, Y=io$Y))
+ # N random starting points gaussian (TODO: around true β?)
+ res <- matrix(nrow=d*K+1, ncol=N)
+ for (i in seq_len(N))
+ {
+ β_init <- rnorm(d*K)
+ par <- op$run( c(rep(1/K,K-1), β_init, rep(0,K)) )
+ par <- op$linArgs(par)
+ Qn <- op$f(par)
+ res[,i] = c(Qn, par[K:(K+d*K-1)])
+ }
+ res #TODO: plot this, not just return it...
+}
--- /dev/null
+#' Generate sample inputs-outputs
+#'
+#' Generate input matrix X of size nxd and binary output of size n, where Y is subdivided
+#' into K groups of proportions p. Inside one group, the probability law P(Y=1) is
+#' described by the corresponding column parameter in the matrix β + intercept b.
+#'
+#' @param n Number of individuals
+#' @param p Vector of K-1 populations relative proportions (sum <= 1)
+#' @param β Vectors of model parameters for each population, of size dxK
+#' @param b Vector of intercept values (use rep(0,K) for no intercept)
+#' @param link Link type; "logit" or "probit"
+#'
+#' @return A list with
+#' \itemize{
+#' \item{X: the input matrix (size nxd)}
+#' \item{Y: the output vector (size n)}
+#' \item{index: the population index (in 1:K) for each row in X}
+#' }
+#'
+#' @export
+generateSampleIO = function(n, p, β, b, link)
+{
+ # Check arguments
+ tryCatch({n = as.integer(n)}, error=function(e) stop("Cannot convert n to integer"))
+ if (length(n) > 1)
+ warning("n is a vector but should be scalar: only first element used")
+ if (n <= 0)
+ stop("n: positive integer")
+ if (!is.matrix(β) || !is.numeric(β) || any(is.na(β)))
+ stop("β: real matrix, no NAs")
+ K = ncol(β)
+ if (!is.numeric(p) || length(p)!=K-1 || any(is.na(p)) || any(p<0) || sum(p) > 1)
+ stop("p: positive vector of size K-1, no NA, sum<=1")
+ p <- c(p, 1-sum(p))
+ if (!is.numeric(b) || length(b)!=K || any(is.na(b)))
+ stop("b: real vector of size K, no NA")
+
+ #random generation of the size of each population in X~Y (unordered)
+ classes = rmultinom(1, n, p)
+
+ d = nrow(β)
+ zero_mean = rep(0,d)
+ id_sigma = diag(rep(1,d))
+ # Always consider an intercept (use b=0 for none)
+ d = d + 1
+ β = rbind(β, b)
+ X = matrix(nrow=0, ncol=d)
+ Y = c()
+ index = c()
+ for (i in 1:ncol(β))
+ {
+ index = c(index, rep(i, classes[i]))
+ newXblock = cbind( MASS::mvrnorm(classes[i], zero_mean, id_sigma), 1 )
+ arg_link = newXblock%*%β[,i]
+ probas =
+ if (link == "logit")
+ {
+ e_arg_link = exp(arg_link)
+ e_arg_link / (1 + e_arg_link)
+ }
+ else #"probit"
+ pnorm(arg_link)
+ probas[is.nan(probas)] = 1 #overflow of exp(x)
+ X = rbind(X, newXblock)
+ Y = c( Y, vapply(probas, function(p) (rbinom(1,1,p)), 1) )
+ }
+ shuffle = sample(n)
+ # Returned X should not contain an intercept column (it's an argument of estimation
+ # methods)
+ list("X"=X[shuffle,-d], "Y"=Y[shuffle], "index"=index[shuffle])
+}
--- /dev/null
+#' normalize
+#'
+#' Normalize a vector or a matrix (by columns), using euclidian norm
+#'
+#' @param X Vector or matrix to be normalized
+#'
+#' @return The normalized matrix (1 column if X is a vector)
+#'
+#' @export
+normalize = function(X)
+{
+ X = as.matrix(X)
+ norm2 = sqrt( colSums(X^2) )
+ sweep(X, 2, norm2, '/')
+}
+
+# Computes a tensor-vector product
+#
+# @param Te third-order tensor (size dxdxd)
+# @param w vector of size d
+#
+# @return Matrix of size dxd
+#
+.T_I_I_w = function(Te, w)
+{
+ d = length(w)
+ Ma = matrix(0,nrow=d,ncol=d)
+ for (j in 1:d)
+ Ma = Ma + w[j] * Te[,,j]
+ Ma
+}
+
+# Computes the second-order empirical moment between input X and output Y
+#
+# @param X matrix of covariates (of size n*d)
+# @param Y vector of responses (of size n)
+#
+# @return Matrix of size dxd
+#
+.Moments_M2 = function(X, Y)
+{
+ n = nrow(X)
+ d = ncol(X)
+ M2 = matrix(0,nrow=d,ncol=d)
+ matrix( .C("Moments_M2", X=as.double(X), Y=as.double(Y), pn=as.integer(n),
+ pd=as.integer(d), M2=as.double(M2), PACKAGE="morpheus")$M2, nrow=d, ncol=d)
+}
+
+# Computes the third-order empirical moment between input X and output Y
+#
+# @param X matrix of covariates (of size n*d)
+# @param Y vector of responses (of size n)
+#
+# @return Array of size dxdxd
+#
+.Moments_M3 = function(X, Y)
+{
+ n = nrow(X)
+ d = ncol(X)
+ M3 = array(0,dim=c(d,d,d))
+ array( .C("Moments_M3", X=as.double(X), Y=as.double(Y), pn=as.integer(n),
+ pd=as.integer(d), M3=as.double(M3), PACKAGE="morpheus")$M3, dim=c(d,d,d) )
+}
+
+#' computeMoments
+#'
+#' Compute cross-moments of order 1,2,3 from X,Y
+#'
+#' @inheritParams computeMu
+#'
+#' @return A list L where L[[i]] is the i-th cross-moment
+#'
+#' @export
+computeMoments = function(X, Y)
+ list( colMeans(Y * X), .Moments_M2(X,Y), .Moments_M3(X,Y) )
+
+# Find the optimal assignment (permutation) between two sets (minimize cost)
+#
+# @param distances The distances matrix, in columns (distances[i,j] is distance between i
+# and j)
+#
+# @return A permutation minimizing cost
+#
+.hungarianAlgorithm = function(distances)
+{
+ n = nrow(distances)
+ .C("hungarianAlgorithm", distances=as.double(distances), pn=as.integer(n),
+ assignment=integer(n), PACKAGE="morpheus")$assignment
+}
+
+#' alignMatrices
+#'
+#' Align a set of parameters matrices, with potential permutations.
+#'
+#' @param Ms A list of matrices, all of same size DxK
+#' @param ref Either a reference matrix or "mean" to align on empirical mean
+#' @param ls_mode How to compute the labels assignment: "exact" for exact algorithm
+#' (default, but might be time-consuming, complexity is O(K^3) ), or "approx1", or
+#' "approx2" to apply a greedy matching algorithm (heuristic) which for each column in
+#' reference (resp. in current row) compare to all unassigned columns in current row
+#' (resp. in reference)
+#'
+#' @return The aligned list (of matrices), of same size as Ms
+#'
+#' @export
+alignMatrices = function(Ms, ref, ls_mode)
+{
+ if (!is.matrix(ref) && ref != "mean")
+ stop("ref: matrix or 'mean'")
+ if (!ls_mode %in% c("exact","approx1","approx2"))
+ stop("ls_mode in {'exact','approx1','approx2'}")
+
+ K <- ncol(Ms[[1]])
+ if (is.character(ref)) #ref=="mean"
+ m_sum = Ms[[1]]
+ L <- length(Ms)
+ for (i in ifelse(is.character(ref),2,1):L)
+ {
+ m_ref = if (is.character(ref)) m_sum / (i-1) else ref
+ m = Ms[[i]] #shorthand
+
+ if (ls_mode == "exact")
+ {
+ #distances[i,j] = distance between m column i and ref column j
+ distances = apply( m_ref, 2, function(col) ( sqrt(colSums((m-col)^2)) ) )
+ assignment = .hungarianAlgorithm(distances)
+ col <- m[,assignment]
+ if (is.list(Ms)) Ms[[i]] <- col else Ms[,,i] <- col
+ }
+ else
+ {
+ # Greedy matching:
+ # approx1: li[[i]][,j] is assigned to m[,k] minimizing dist(li[[i]][,j],m[,k'])
+ # approx2: m[,j] is assigned to li[[i]][,k] minimizing dist(m[,j],li[[i]][,k'])
+ available_indices = 1:K
+ for (j in 1:K)
+ {
+ distances =
+ if (ls_mode == "approx1")
+ {
+ apply(as.matrix(m[,available_indices]), 2,
+ function(col) ( sqrt(sum((col - m_ref[,j])^2)) ) )
+ }
+ else #approx2
+ {
+ apply(as.matrix(m_ref[,available_indices]), 2,
+ function(col) ( sqrt(sum((col - m[,j])^2)) ) )
+ }
+ indMin = which.min(distances)
+ if (ls_mode == "approx1")
+ {
+ col <- m[ , available_indices[indMin] ]
+ if (is.list(Ms)) Ms[[i]][,j] <- col else Ms[,j,i] <- col
+ }
+ else #approx2
+ {
+ col <- available_indices[indMin]
+ if (is.list(Ms)) Ms[[i]][,col] <- m[,j] else Ms[,col,i] <- m[,j]
+ }
+ available_indices = available_indices[-indMin]
+ }
+ }
+
+ # Update current sum with "label-switched" li[[i]]
+ if (is.character(ref)) #ref=="mean"
+ m_sum = m_sum + Ms[[i]]
+ }
+ Ms
+}
--- /dev/null
+\name{morpheus-package}
+\alias{morpheus-package}
+\alias{morpheus}
+\docType{package}
+
+\title{
+ \packageTitle{morpheus}
+}
+
+\description{
+ \packageDescription{morpheus}
+}
+
+\details{
+ The package devtools should be useful in development stage, since we rely on testthat for
+ unit tests, and roxygen2 for documentation. knitr is used to generate the package vignette.
+ Concerning the other suggested packages:
+ \itemize{
+ \item{tensor is used for comparing to some reference functions initially coded in R;
+ it should not be required in further package versions;}
+ \item{jointDiag allows to solve a joint diagonalization problem, providing a more
+ robust solution compared to a single diagonalization;}
+ \item{parallel (generally) permits to run the bootstrap method faster.}
+ }
+
+ The three main functions are located in R/main.R:
+ \itemize{
+ \item{getParamsDirs_ref: reference method to estimate parameters directions;}
+ \item{getParamsDirs: method of choice to estimate parameters directions, using a
+ spectral decomposition of inputs/outputs;}
+ \item{getBootstrapParams: run getParamsDirs on B bootstrap replicates.}
+ }
+}
+
+\author{
+ \packageAuthor{morpheus}
+
+ Maintainer: \packageMaintainer{morpheus}
+}
+
+%\references{
+% TODO: Literature or other references for background information
+%}
+
+%\examples{
+% TODO: simple examples of the most important functions
+%}
--- /dev/null
+#include <stdlib.h>
+
+//index matrix (by columns)
+int mi(int i, int j, int d1, int d2)
+{
+ return j*d1+i;
+}
+
+//index 3-tensor (by columns, matrices ordered by last dim)
+int ti(int i, int j, int k, int d1, int d2, int d3)
+{
+ return k*d1*d2 + j*d1 + i;
+}
+
+// Emprical cross-moment of order 2 between X size nxd and Y size n
+void Moments_M2(double* X, double* Y, int* pn, int* pd, double* M2)
+{
+ int n=*pn, d=*pd;
+ //double* M2 = (double*)calloc(d*d,sizeof(double));
+
+ // M2 = E[Y*X^*2] - E[Y*e^*2] = E[Y (X^*2 - I)]
+ for (int j=0; j<d; j++)
+ {
+ for (int i=0; i<n; i++)
+ {
+ M2[mi(j,j,d,d)] -= Y[i] / n;
+ for (int k=0; k<d; k++)
+ M2[mi(j,k,d,d)] += Y[i] * X[mi(i,j,n,d)]*X[mi(i,k,n,d)] / n;
+ }
+ }
+}
+
+// Emprical cross-moment of order 3 between X size nxd and Y size n
+void Moments_M3(double* X, double* Y, int* pn, int* pd, double* M3)
+{
+ int n=*pn, d=*pd;
+ //double* M3 = (double*)calloc(d*d*d,sizeof(double));
+
+ // M3 = E[Y*X^*3] - E[Y*e*X*e] - E[Y*e*e*X] - E[Y*X*e*e]
+ for (int j=0; j<d; j++)
+ {
+ for (int k=0; k<d; k++)
+ {
+ for (int i=0; i<n; i++)
+ {
+ double tensor_elt = Y[i]*X[mi(i,k,n,d)] / n;
+ M3[ti(j,k,j,d,d,d)] -= tensor_elt;
+ M3[ti(j,j,k,d,d,d)] -= tensor_elt;
+ M3[ti(k,j,j,d,d,d)] -= tensor_elt;
+ for (int o=0; o<d; o++)
+ M3[ti(j,k,o,d,d,d)] += Y[i] * X[mi(i,j,n,d)]*X[mi(i,k,n,d)]*X[mi(i,o,n,d)] / n;
+ }
+ }
+ }
+}
--- /dev/null
+/********************************************************************
+ ********************************************************************
+ **
+ ** libhungarian by Cyrill Stachniss, 2004
+ ** Modified by Benjamin Auder to work on floating point number,
+ ** and to fit into a R package (2017)
+ **
+ ** Solving the Minimum Assignment Problem using the Hungarian Method.
+ **
+ ** ** This file may be freely copied and distributed! **
+ **
+ ** Parts of the used code was originally provided by the
+ ** "Stanford GraphGase", but I made changes to this code.
+ ** As asked by the copyright node of the "Stanford GraphGase",
+ ** I hereby proclaim that this file are *NOT* part of the
+ ** "Stanford GraphGase" distrubition!
+ **
+ ** This file is distributed in the hope that it will be useful,
+ ** but WITHOUT ANY WARRANTY; without even the implied
+ ** warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
+ ** PURPOSE.
+ **
+ ********************************************************************
+ ********************************************************************/
+
+#include <stdio.h>
+#include <stdlib.h>
+
+#define HUNGARIAN_NOT_ASSIGNED 0
+#define HUNGARIAN_ASSIGNED 1
+
+#define HUNGARIAN_MODE_MINIMIZE_COST 0
+#define HUNGARIAN_MODE_MAXIMIZE_UTIL 1
+
+#define INF (1e32)
+#define verbose (0)
+
+typedef struct {
+ int num_rows;
+ int num_cols;
+ double** cost;
+ int** assignment;
+} hungarian_problem_t;
+
+int hungarian_imax(int a, int b) {
+ return (a<b)?b:a;
+}
+
+/** This method initialize the hungarian_problem structure and init
+ * the cost matrices (missing lines or columns are filled with 0).
+ * It returns the size of the quadratic(!) assignment matrix. **/
+int hungarian_init(hungarian_problem_t* p, double** cost_matrix, int rows, int cols,
+ int mode)
+{
+ int i,j;
+ double max_cost = 0.;
+ int org_cols = cols,
+ org_rows = rows;
+
+ // is the number of cols not equal to number of rows ?
+ // if yes, expand with 0-cols / 0-cols
+ rows = hungarian_imax(cols, rows);
+ cols = rows;
+
+ p->num_rows = rows;
+ p->num_cols = cols;
+
+ p->cost = (double**)calloc(rows,sizeof(double*));
+ p->assignment = (int**)calloc(rows,sizeof(int*));
+
+ for(i=0; i<p->num_rows; i++) {
+ p->cost[i] = (double*)calloc(cols,sizeof(double));
+ p->assignment[i] = (int*)calloc(cols,sizeof(int));
+ for(j=0; j<p->num_cols; j++) {
+ p->cost[i][j] = (i < org_rows && j < org_cols) ? cost_matrix[i][j] : 0.;
+ p->assignment[i][j] = 0;
+
+ if (max_cost < p->cost[i][j])
+ max_cost = p->cost[i][j];
+ }
+ }
+
+ if (mode == HUNGARIAN_MODE_MAXIMIZE_UTIL) {
+ for(i=0; i<p->num_rows; i++) {
+ for(j=0; j<p->num_cols; j++)
+ p->cost[i][j] = max_cost - p->cost[i][j];
+ }
+ }
+ else if (mode == HUNGARIAN_MODE_MINIMIZE_COST) {
+ // nothing to do
+ }
+// else
+// fprintf(stderr,"%s: unknown mode. Mode was set to \
+// HUNGARIAN_MODE_MINIMIZE_COST !\n", __FUNCTION__);
+
+ return rows;
+}
+
+/** Free the memory allocated by init. **/
+void hungarian_free(hungarian_problem_t* p) {
+ int i;
+ for(i=0; i<p->num_rows; i++) {
+ free(p->cost[i]);
+ free(p->assignment[i]);
+ }
+ free(p->cost);
+ free(p->assignment);
+ p->cost = NULL;
+ p->assignment = NULL;
+}
+
+/** This method computes the optimal assignment. **/
+void hungarian_solve(hungarian_problem_t* p)
+{
+ int i, j, m, n, k, l, t, q, unmatched;
+ double cost, s;
+ int* col_mate;
+ int* row_mate;
+ int* parent_row;
+ int* unchosen_row;
+ double* row_dec;
+ double* col_inc;
+ double* slack;
+ int* slack_row;
+
+ cost = 0.;
+ m =p->num_rows;
+ n =p->num_cols;
+
+ col_mate = (int*)calloc(p->num_rows,sizeof(int));
+ unchosen_row = (int*)calloc(p->num_rows,sizeof(int));
+ row_dec = (double*)calloc(p->num_rows,sizeof(double));
+ slack_row = (int*)calloc(p->num_rows,sizeof(int));
+
+ row_mate = (int*)calloc(p->num_cols,sizeof(int));
+ parent_row = (int*)calloc(p->num_cols,sizeof(int));
+ col_inc = (double*)calloc(p->num_cols,sizeof(double));
+ slack = (double*)calloc(p->num_cols,sizeof(double));
+
+ for (i=0; i<p->num_rows; i++) {
+ col_mate[i]=0;
+ unchosen_row[i]=0;
+ row_dec[i]=0.;
+ slack_row[i]=0;
+ }
+ for (j=0; j<p->num_cols; j++) {
+ row_mate[j]=0;
+ parent_row[j] = 0;
+ col_inc[j]=0.;
+ slack[j]=0.;
+ }
+
+ for (i=0; i<p->num_rows; ++i)
+ for (j=0; j<p->num_cols; ++j)
+ p->assignment[i][j]=HUNGARIAN_NOT_ASSIGNED;
+
+ // Begin subtract column minima in order to start with lots of zeroes 12
+ for (l=0; l<n; l++)
+ {
+ s=p->cost[0][l];
+ for (k=1; k<m; k++)
+ if (p->cost[k][l]<s)
+ s=p->cost[k][l];
+ cost+=s;
+ if (s!=0.)
+ for (k=0; k<m; k++)
+ p->cost[k][l]-=s;
+ }
+ // End subtract column minima in order to start with lots of zeroes 12
+
+ // Begin initial state 16
+ t=0;
+ for (l=0; l<n; l++)
+ {
+ row_mate[l]= -1;
+ parent_row[l]= -1;
+ col_inc[l]=0.;
+ slack[l]=INF;
+ }
+ for (k=0; k<m; k++)
+ {
+ s=p->cost[k][0];
+ for (l=1; l<n; l++)
+ if (p->cost[k][l]<s)
+ s=p->cost[k][l];
+ row_dec[k]=s;
+ for (l=0; l<n; l++)
+ if (s==p->cost[k][l] && row_mate[l]<0)
+ {
+ col_mate[k]=l;
+ row_mate[l]=k;
+ goto row_done;
+ }
+ col_mate[k]= -1;
+ unchosen_row[t++]=k;
+row_done:
+ ;
+ }
+ // End initial state 16
+
+ // Begin Hungarian algorithm 18
+ if (t==0)
+ goto done;
+ unmatched=t;
+ while (1)
+ {
+ q=0;
+ while (1)
+ {
+ while (q<t)
+ {
+ // Begin explore node q of the forest 19
+ {
+ k=unchosen_row[q];
+ s=row_dec[k];
+ for (l=0; l<n; l++)
+ if (slack[l] > 0.)
+ {
+ double del=p->cost[k][l]-s+col_inc[l];
+ if (del<slack[l])
+ {
+ if (del==0.)
+ {
+ if (row_mate[l]<0)
+ goto breakthru;
+ slack[l]=0.;
+ parent_row[l]=k;
+ unchosen_row[t++]=row_mate[l];
+ }
+ else
+ {
+ slack[l]=del;
+ slack_row[l]=k;
+ }
+ }
+ }
+ }
+ // End explore node q of the forest 19
+ q++;
+ }
+
+ // Begin introduce a new zero into the matrix 21
+ s=INF;
+ for (l=0; l<n; l++)
+ if (slack[l]>0. && slack[l]<s)
+ s=slack[l];
+ for (q=0; q<t; q++)
+ row_dec[unchosen_row[q]]+=s;
+ for (l=0; l<n; l++)
+ if (slack[l]>0.)
+ {
+ slack[l]-=s;
+ if (slack[l]==0.)
+ {
+ // Begin look at a new zero 22
+ k=slack_row[l];
+ if (row_mate[l]<0)
+ {
+ for (j=l+1; j<n; j++)
+ if (slack[j]==0.)
+ col_inc[j]+=s;
+ goto breakthru;
+ }
+ else
+ {
+ parent_row[l]=k;
+ unchosen_row[t++]=row_mate[l];
+ }
+ // End look at a new zero 22
+ }
+ }
+ else
+ col_inc[l]+=s;
+ // End introduce a new zero into the matrix 21
+ }
+breakthru:
+ // Begin update the matching 20
+ while (1)
+ {
+ j=col_mate[k];
+ col_mate[k]=l;
+ row_mate[l]=k;
+ if (j<0)
+ break;
+ k=parent_row[j];
+ l=j;
+ }
+ // End update the matching 20
+ if (--unmatched==0)
+ goto done;
+ // Begin get ready for another stage 17
+ t=0;
+ for (l=0; l<n; l++)
+ {
+ parent_row[l]= -1;
+ slack[l]=INF;
+ }
+ for (k=0; k<m; k++)
+ if (col_mate[k]<0)
+ unchosen_row[t++]=k;
+ // End get ready for another stage 17
+ }
+done:
+
+/* // Begin doublecheck the solution 23
+ for (k=0; k<m; k++)
+ for (l=0; l<n; l++)
+ if (p->cost[k][l]<row_dec[k]-col_inc[l])
+ exit(0);
+ for (k=0; k<m; k++)
+ {
+ l=col_mate[k];
+ if (l<0 || p->cost[k][l]!=row_dec[k]-col_inc[l])
+ exit(0);
+ }
+ k=0;
+ for (l=0; l<n; l++)
+ if (col_inc[l])
+ k++;
+ if (k>m)
+ exit(0);
+ // End doublecheck the solution 23
+*/ // End Hungarian algorithm 18
+
+ for (i=0; i<m; ++i)
+ {
+ p->assignment[i][col_mate[i]]=HUNGARIAN_ASSIGNED;
+ /*TRACE("%d - %d\n", i, col_mate[i]);*/
+ }
+ for (k=0; k<m; ++k)
+ {
+ for (l=0; l<n; ++l)
+ {
+ /*TRACE("%d ",p->cost[k][l]-row_dec[k]+col_inc[l]);*/
+ p->cost[k][l]=p->cost[k][l]-row_dec[k]+col_inc[l];
+ }
+ /*TRACE("\n");*/
+ }
+ for (i=0; i<m; i++)
+ cost+=row_dec[i];
+ for (i=0; i<n; i++)
+ cost-=col_inc[i];
+
+ free(slack);
+ free(col_inc);
+ free(parent_row);
+ free(row_mate);
+ free(slack_row);
+ free(row_dec);
+ free(unchosen_row);
+ free(col_mate);
+}
+
+
+/***************
+ * Usage in R :
+ **************/
+
+// Auxiliary function for hungarianAlgorithm below
+double** array_to_matrix(double* m, int rows, int cols)
+{
+ double** r = (double**)calloc(rows,sizeof(double*));
+ for (int i=0; i<rows; i++)
+ {
+ r[i] = (double*)calloc(cols,sizeof(double));
+ for (int j=0; j<cols; j++)
+ r[i][j] = m[i*cols+j];
+ }
+ return r;
+}
+
+//TODO: re-code this algorithm in a more readable way, based on
+//https://www.topcoder.com/community/data-science/data-science-tutorials/\
+// assignment-problem-and-hungarian-algorithm/
+// Get the optimal assignment, by calling hungarian_solve above; "distances" in columns
+void hungarianAlgorithm(double* distances, int* pn, int* assignment)
+{
+ int n = *pn;
+ double** mat = array_to_matrix(distances, n, n);
+
+ hungarian_problem_t p;
+ hungarian_init(&p, mat, n, n, HUNGARIAN_MODE_MINIMIZE_COST);
+
+ hungarian_solve(&p);
+
+ // Copy the optimal assignment in a R-friendly structure
+ for (int i=0; i < n; i++) {
+ for (int j=0; j < n; j++) {
+ if (p.assignment[i][j])
+ assignment[i] = j+1; //add 1 since we work with R
+ }
+ }
+
+ hungarian_free(&p);
+ for (int i=0; i<n; i++)
+ free(mat[i]);
+ free(mat);
+}
--- /dev/null
+library(testthat)
+#library(morpheus)
+load_all()
+
+test_check("morpheus")
--- /dev/null
+context("Moments_M2")
+
+# (Slower, but trusted) R version of Moments_M2
+.Moments_M2_R = function(X,Y)
+{
+ library(tensor)
+ n = nrow(X)
+ d = ncol(X)
+ v = rep(0,d)
+ e = diag(rep(1,d))
+
+ M21 = M2_1 = tensor(v,v)
+ for (i in 1:n)
+ M21 = M21 + Y[i] * tensor(X[i,],X[i,])
+ M21 = (1/n) * M21
+
+ for (j in 1:d)
+ {
+ L = tensor(v,v)
+ for (i in 1:n)
+ L = L + Y[i]*tensor(e[,j],e[,j])
+ L = (1/n) * L
+ M2_1 = M2_1 + L
+ }
+
+ M2 = M21 - M2_1
+ return (M2)
+}
+
+test_that("both versions of Moments_M2 agree on various inputs",
+{
+ for (n in c(20,200))
+ {
+ for (d in c(2,10,20))
+ {
+ X = matrix( runif(n*d,min=-1,max=1), nrow=n )
+ Y = runif(n,min=-1,max=1)
+ M2 = .Moments_M2(X,Y)
+ M2_R = .Moments_M2_R(X,Y)
+ expect_equal(max(abs(M2 - M2_R)), 0)
+ }
+ }
+})
--- /dev/null
+context("Moments_M3")
+
+# (Slower, but trusted) R version of Moments_M3
+.Moments_M3_R = function(X,Y)
+{
+ library(tensor)
+ n = nrow(X)
+ d = ncol(X)
+ v = rep(0,d)
+ e = diag(rep(1,d))
+
+ M31=M3_1=M3_2=M3_3 = tensor(tensor(v,v),v)
+ for (i in 1:n)
+ M31 = M31 + (Y[i]*tensor(tensor(X[i,],X[i,]),X[i,]))
+ M31 = (1/n)*M31
+
+ for(j in 1:d)
+ {
+ l1=l2=l3 = tensor(tensor(v,v),v)
+ for(i in 1:n)
+ {
+ l1 = l1 + Y[i]*tensor(tensor(e[,j],X[i,]),e[,j])
+ l2 = l2 + Y[i]*tensor(tensor(e[,j],e[,j]),X[i,])
+ l3 = l3 + Y[i]*tensor(tensor(X[i,],e[,j]),e[,j])
+ }
+ l1 = (1/n)*l1
+ l2 = (1/n)*l2
+ l3 = (1/n)*l3
+ M3_1 = M3_1+l1
+ M3_2 = M3_2+l2
+ M3_3 = M3_3+l3
+ }
+
+ M3 = M31-(M3_1+M3_2+M3_3)
+ return (M3)
+}
+
+test_that("both versions of Moments_M3 agree on various inputs",
+{
+ for (n in c(20,200))
+ {
+ for (d in c(2,10,20))
+ {
+ X = matrix( runif(n*d,min=-1,max=1), nrow=n )
+ Y = runif(n,min=-1,max=1)
+ M3 = .Moments_M3(X,Y)
+ M3_R = .Moments_M3_R(X,Y)
+ expect_equal(max(abs(M3 - M3_R)), 0)
+ }
+ }
+})
--- /dev/null
+context("alignMatrices")
+
+# Helper to generate a random series of matrices to align
+.generateMatrices = function(d, K, N, noise)
+{
+ matrices = list( matrix(runif(d*K, min=-1, max=1),ncol=K) ) #reference
+ for (i in 2:N)
+ {
+ matrices[[i]] <- matrices[[1]][,sample(1:K)]
+ if (noise)
+ matrices[[i]] = matrices[[i]] + matrix(rnorm(d*K, sd=0.05), ncol=K)
+ }
+ matrices
+}
+
+test_that("labelSwitchingAlign correctly aligns de-noised parameters",
+{
+ N = 30 #number of matrices
+ d_K_list = list(c(2,2), c(5,3))
+ for (i in 1:2)
+ {
+ d = d_K_list[[i]][1]
+ K = d_K_list[[i]][2]
+
+ # 1] Generate matrix series
+ matrices_permut = .generateMatrices(d,K,N,noise=FALSE)
+
+ # 2] Call align function with mode=approx1
+ matrices_aligned =
+ alignMatrices(matrices_permut[2:N], ref=matrices_permut[[1]], ls_mode="approx1")
+
+ # 3] Check alignment
+ for (j in 2:N)
+ expect_equal(matrices_aligned[[j-1]], matrices_permut[[1]])
+
+ # 2bis] Call align function with mode=approx2
+ matrices_aligned =
+ alignMatrices(matrices_permut[2:N], ref=matrices_permut[[1]], ls_mode="approx2")
+
+ # 3bis] Check alignment
+ for (j in 2:N)
+ expect_equal(matrices_aligned[[j-1]], matrices_permut[[1]])
+ }
+})
+
+test_that("labelSwitchingAlign correctly aligns noisy parameters",
+{
+ N = 30 #number of matrices
+ d_K_list = list(c(2,2), c(5,3))
+ for (i in 1:2)
+ {
+ d = d_K_list[[i]][1]
+ K = d_K_list[[i]][2]
+ max_error = d * 0.2 #TODO: what value to choose ?
+
+ # 1] Generate matrix series
+ matrices_permut = .generateMatrices(d,K,N,noise=TRUE)
+
+ # 2] Call align function
+ matrices_aligned = alignMatrices(matrices_permut, ref="mean", ls_mode="exact")
+
+ # 3] Check alignment
+ for (j in 2:N)
+ {
+ expect_that( norm(matrices_aligned[[j]] - matrices_permut[[1]]),
+ is_less_than(max_error) )
+ }
+ }
+})
--- /dev/null
+context("computeMu")
+
+test_that("on input of sufficient size, β/||β|| is estimated accurately enough",
+{
+ n = 100000
+ d = 2
+ K = 2
+ p = 1/2
+
+ βs_ref = array( c(1,0,0,1 , 1,-2,3,1), dim=c(d,K,2) )
+ for (i in 1:(dim(βs_ref)[3]))
+ {
+ μ_ref = normalize(βs_ref[,,i])
+ for (model in c("logit","probit"))
+ {
+ cat("\n\n",model," :\n",sep="")
+
+ io = generateSampleIO(n, p, βs_ref[,,i], rep(0,K), model)
+ μ = computeMu(io$X, io$Y, list(K=K))
+ μ_aligned = alignMatrices(list(μ), ref=μ_ref, ls_mode="exact")[[1]]
+
+ #Some traces: 0 is not well estimated, but others are OK
+ cat("Reference normalized matrix:\n")
+ print(μ_ref)
+ cat("Estimated normalized matrix:\n")
+ print(μ_aligned)
+ cat("Difference norm (Matrix norm ||.||_1, max. abs. sum on a column)\n")
+ diff_norm = norm(μ_ref - μ_aligned)
+ cat(diff_norm,"\n")
+
+ #NOTE: 0.5 is loose threshold, but values around 0.3 are expected...
+ expect_that( diff_norm, is_less_than(0.5) )
+ }
+ }
+})
--- /dev/null
+context("hungarianAlgorithm")
+
+test_that("HungarianAlgorithm provides the correct answer on various inputs",
+{
+ for (n in c(2,3,10,50))
+ {
+ for (repetition in 1:3)
+ {
+ # Generate a random vector, and permute it:
+ # we expect the algorithm to retrieve the permutation
+ v = runif(n, min=-1, max=1)
+ permutation = sample(1:n)
+ v_p = v[permutation]
+
+ # v is reference, v_p is v after permutation
+ # distances[i,j] = distance between i-th component of v and j-th component of v_p
+ # "in rows : reference; in columns: after permutation"
+ # "permutation" contains order on v to match v_p as closely as possible
+ distances = sapply(v_p, function(elt) ( abs(elt - v) ) )
+ assignment = .hungarianAlgorithm(distances)
+ expect_equal(assignment, permutation)
+ }
+ }
+})
--- /dev/null
+context("jointDiag::ajd")
+
+#auxiliary to test diagonality
+.computeMuCheckDiag = function(X, Y, K, jd_method, β_ref)
+{
+ d = ncol(X)
+ #TODO: redundant code, same as computeMu() main method. Comments are stripped
+ M3 = .Moments_M3(X,Y)
+ M2_t = array(dim=c(d,d,d))
+ for (i in 1:d)
+ {
+ ρ = c(rep(0,i-1),1,rep(0,d-i))
+ M2_t[,,i] = .T_I_I_w(M3,ρ)
+ }
+ jd = jointDiag::ajd(M2_t, method=jd_method)
+ V = if (jd_method=="uwedge") jd$B else solve(jd$A)
+ M2_t = array(dim=c(d,d,K))
+ for (i in 1:K)
+ M2_t[,,i] = .T_I_I_w(M3,V[,i])
+ #END of computeMu() code
+
+ max_error = 0.5 #TODO: tune ?
+ invβ = MASS::ginv(β_ref)
+ for (i in 1:K)
+ {
+ shouldBeDiag = invβ %*% M2_t[,,i] %*% t(invβ)
+ expect_that(
+ mean( abs(shouldBeDiag[upper.tri(shouldBeDiag) | lower.tri(shouldBeDiag)]) ),
+ is_less_than(max_error) )
+ }
+}
+
+test_that("'jedi' and 'uwedge' joint-diagonalization methods return a correct matrix",
+{
+ n = 10000
+ d_K = list( c(2,2), c(5,3), c(20,13) )
+
+ for (dk_index in 1:length(d_K))
+ {
+ d = d_K[[dk_index]][1]
+ K = d_K[[dk_index]][2]
+ #NOTE: sometimes large errors if pr is not balanced enough (e.g. random);
+ # same note for β. However we could be more random than that...
+ β_ref = rbind(diag(K),matrix(0,nrow=d-K,ncol=K))
+ io = generateSampleIO(n, p=rep(1/K,K-1), β=β_ref, rep(0,K), link="logit")
+ .computeMuCheckDiag(io$X, io$Y, K, jd_method="uwedge", β_ref)
+ #TODO: some issues with jedi method (singular system)
+ #.computeMuCheckDiag(io$X, io$Y, K, jd_method="jedi", β_ref)
+ }
+})
--- /dev/null
+context("OptimParams")
+
+naive_f = function(link, M1,M2,M3, p,β,b)
+{
+ d = length(M1)
+ K = length(p)
+ λ <- sqrt(colSums(β^2))
+
+ # Compute β x2,3 (self) tensorial products
+ β2 = array(0, dim=c(d,d,K))
+ β3 = array(0, dim=c(d,d,d,K))
+ for (k in 1:K)
+ {
+ for (i in 1:d)
+ {
+ for (j in 1:d)
+ {
+ β2[i,j,k] = β[i,k]*β[j,k]
+ for (l in 1:d)
+ β3[i,j,l,k] = β[i,k]*β[j,k]*β[l,k]
+ }
+ }
+ }
+
+ res = 0
+ for (i in 1:d)
+ {
+ term = 0
+ for (k in 1:K)
+ term = term + p[k]*.G(link,1,λ[k],b[k])*β[i,k]
+ res = res + (term - M1[i])^2
+ for (j in 1:d)
+ {
+ term = 0
+ for (k in 1:K)
+ term = term + p[k]*.G(link,2,λ[k],b[k])*β2[i,j,k]
+ res = res + (term - M2[i,j])^2
+ for (l in 1:d)
+ {
+ term = 0
+ for (k in 1:K)
+ term = term + p[k]*.G(link,3,λ[k],b[k])*β3[i,j,l,k]
+ res = res + (term - M3[i,j,l])^2
+ }
+ }
+ }
+ res
+}
+
+test_that("naive computation provides the same result as vectorized computations",
+{
+ h <- 1e-7 #for finite-difference tests
+ tol <- .25 * sqrt(h) #about 7.9 e-5
+ for (dK in list( c(2,2), c(5,3)))
+ {
+ d = dK[1]
+ K = dK[2]
+
+ M1 = runif(d, -1, 1)
+ M2 = matrix(runif(d*d,-1,1), ncol=d)
+ M3 = array(runif(d*d*d,-1,1), dim=c(d,d,d))
+
+ for (link in c("logit","probit"))
+ {
+ op = new("OptimParams", "li"=link, "M1"=as.double(M1),
+ "M2"=as.double(M2), "M3"=as.double(M3), "K"=as.integer(K))
+
+ for (var in seq_len((2+d)*K-1))
+ {
+ p = runif(K, 0, 1)
+ p = p / sum(p)
+ β <- matrix(runif(d*K,-5,5),ncol=K)
+ b = runif(K, -5, 5)
+ x <- c(p[1:(K-1)],as.double(β),b)
+
+ # Test functions values
+ expect_equal( op$f(x), naive_f(link,M1,M2,M3, p,β,b) )
+
+ # Test finite differences ~= gradient values
+ dir_h <- rep(0, (2+d)*K-1)
+ dir_h[var] = h
+
+ expect_equal( op$grad_f(x)[var], ( op$f(x+dir_h) - op$f(x) ) / h, tol )
+ }
+ }
+ }
+})
--- /dev/null
+#!/bin/bash
+
+# Synchronize package folder
+rsync -a --delete pkg/ pkg-cran/
+cd pkg-cran
+
+# Replace all used greek letters by beta, mu, ...etc
+for file in `find -name *.R`; do
+ if [ -f $file ]; then
+ sed -i 's/μ/mu/g' $file
+ sed -i 's/β/beta/g' $file
+ sed -i 's/λ/lambda/g' $file
+ sed -i 's/Σ/Sigma/g' $file
+ fi
+done
+
+if [ "$1" == "i" ]; then
+ # Install package
+ echo "rp()" | R --slave
+fi
+
+cd ..
projects / morpheus.git / commitdiff
