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.
+Warning: tests won't run well on the non-UTF8 package. Run ./to-cran.sh and
+then inside pkg-cran call devtools::test().
+
## Example
Install the package, then ?multiRun is a possible starting point.
--- /dev/null
+#$# git-fat 00f581aca4dd370615fa0ea99e730d6dd42fe347 585618
devtools,
flexmix,
parallel,
- testthat,
+ testthat (>= 3.0.0),
roxygen2
License: MIT + file LICENSE
-RoxygenNote: 7.0.2
+RoxygenNote: 7.1.1
URL: https://github.com/yagu0/morpheus
Collate:
'utils.R'
'optimParams.R'
'plot.R'
'sampleIO.R'
+Config/testthat/edition: 3
#' Estimate the normalized columns μ of the β matrix parameter in a mixture of
#' logistic regressions models, with a spectral method described in the package vignette.
#'
+#' @name computeMu
+#'
#' @param X Matrix of input data (size nxd)
#' @param Y Vector of binary outputs (size n)
#' @param optargs List of optional argument:
#' 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
+#' 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())
+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)
+ 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
+ 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
+ Σ <- svd(M[[2]])$d
large_ratio <- ( abs(Σ[-d] / Σ[-1]) > 3 )
K <- if (any(large_ratio)) max(2, which.min(large_ratio)) else d
}
stop("K: integer >= 2, <= 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)
+ 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
+ jd_nvects <- d
+ fixed_design <- TRUE
}
- M2_t = array(dim=c(d,d,jd_nvects))
+ 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)
+ 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 =
+ 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)})
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))
+ 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])
+ 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") MASS::ginv(jd$B) else jd$A
- μ = normalize(U[,1:K])
+ 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]
+ C <- MASS::ginv(μ) %*% M[[1]]
+ μ[,C < 0] <- - μ[,C < 0]
μ
}
#' basis (exact same inputs). It's even possible to compare methods on some
#' deterministic design of experiments.
#'
+#' @name multiRun
+#'
#' @param fargs List of arguments for the estimation functions
#' @param estimParams List of nf function(s) to apply on fargs
#' @param prepareArgs Prepare arguments for the functions inside estimParams
#' library(flexmix)
#' ind <- fargs$ind
#' K <- fargs$K
-#' dat = as.data.frame( cbind(fargs$Y[ind],fargs$X[ind,]) )
-#' out = refit( flexmix( cbind(V1, 1 - V1) ~ 0+., data=dat, k=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) )
#' } ),
#' res[[i]] <- alignMatrices(res[[i]], ref=μ, ls_mode="exact")
#'
#' # Monte-Carlo + optimParams from X,Y, morpheus VS flexmix
-#' res <- multiRun(list(n=1000,p=1/2,β=β,b=c(0,0),link="logit"),list(
+#' res <- multiRun(list(n=1000,p=1/2,β=β,b=c(0,0),link="logit"), list(
#' # morpheus
#' function(fargs) {
#' library(morpheus)
#' library(flexmix)
#' K <- fargs$K
#' dat <- as.data.frame( cbind(fargs$Y,fargs$X) )
-#' out <- refit( flexmix( cbind(V1, 1 - V1) ~ 0+., data=dat, k=K,
+#' out <- refit( flexmix( cbind(V1, 1 - V1) ~ ., data=dat, k=K,
#' model=FLXMRglm(family="binomial") ) )
#' sapply( seq_len(K), function(i)
-#' as.double( out@@components[[1]][[i]][,1] ) )
+#' as.double( out@@components[[1]][[i]][2:(1+ncol(fargs$X)),1] ) )
#' } ),
#' prepareArgs = function(fargs,index) {
#' library(morpheus)
-#' io = generateSampleIO(fargs$n, fargs$p, fargs$β, fargs$b, fargs$link)
-#' fargs$X = io$X
-#' fargs$Y = io$Y
-#' fargs$K = ncol(fargs$β)
-#' fargs$link = fargs$link
-#' fargs$M = computeMoments(io$X,io$Y)
+#' io <- generateSampleIO(fargs$n, fargs$p, fargs$β, fargs$b, fargs$link)
+#' fargs$X <- io$X
+#' fargs$Y <- io$Y
+#' fargs$K <- ncol(fargs$β)
+#' fargs$link <- fargs$link
+#' fargs$M <- computeMoments(io$X,io$Y)
#' fargs
#' }, N=10, ncores=3)
#' for (i in 1:2)
if (ncores > 1)
{
- cl = parallel::makeCluster(ncores, outfile="")
+ cl <- parallel::makeCluster(ncores, outfile="")
parallel::clusterExport(cl, c("fargs","verbose"), environment())
- list_res = parallel::clusterApplyLB(cl, 1:N, estimParamAtIndex)
+ list_res <- parallel::clusterApplyLB(cl, 1:N, estimParamAtIndex)
parallel::stopCluster(cl)
}
else
- list_res = lapply(1:N, estimParamAtIndex)
+ list_res <- lapply(1:N, estimParamAtIndex)
# De-interlace results: output one list per function
nf <- length(estimParams)
+#' optimParams
+#'
#' Wrapper function for OptimParams class
#'
+#' @name optimParams
+#'
#' @param X Data matrix of covariables
#' @param Y Output as a binary vector
#' @param K Number of populations.
#' # Optimize parameters from estimated μ
#' io <- generateSampleIO(100,
#' 1/2, matrix(c(1,-2,3,1),ncol=2), c(0,0), "logit")
-#' μ = computeMu(io$X, io$Y, list(K=2))
+#' μ <- computeMu(io$X, io$Y, list(K=2))
#' o <- optimParams(io$X, io$Y, 2, "logit")
#' \donttest{
#' θ0 <- list(p=1/2, β=μ, b=c(0,0))
{
stop("θ0$p: length K-1, no NA, positive integers, sum to <= 1")
}
- if (is.null(θ0$b))
+ # NOTE: [["b"]] instead of $b because $b would match $beta (in pkg-cran)
+ if (is.null(θ0[["b"]]))
θ0$b = rep(0, K)
else if (!is.numeric(θ0$b) || length(θ0$b) != K || any(is.na(θ0$b)))
stop("θ0$b: length K, no NA")
# (Re)Set W to identity, to allow several run from the same object
W <<- if (is.null(userW)) diag(d+d^2+d^3) else userW
- #NOTE: loopMax = 3 seems to not improve the final results.
+ # NOTE: loopMax = 3 seems to not improve the final results.
loopMax <- ifelse(is.null(userW), 2, 1)
x_init <- linArgs(θ0)
for (loop in 1:loopMax)
#'
#' Plot compared histograms of a single parameter (scalar)
#'
+#' @name plotHist
+#'
#' @param mr Output of multiRun(), list of lists of functions results
#' @param x Row index of the element inside the aggregated parameter
#' @param y Column index of the element inside the aggregated parameter
plotHist <- function(mr, x, y, ...)
{
params <- .extractParam(mr, x, y)
- L = length(params)
+ 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))
args <- list(...)
}
}
+# NOTE: roxygen2 bug, "@inheritParams plotHist" fails in next header:
+
#' plotBox
#'
#' Draw compared boxplots of a single parameter (scalar)
#'
-#' @inheritParams plotHist
+#' @name plotBox
+#'
+#' @param mr Output of multiRun(), list of lists of functions results
+#' @param x Row index of the element inside the aggregated parameter
+#' @param y Column index of the element inside the aggregated parameter
+#' @param ... Additional graphical parameters (xlab, ylab, ...)
#'
#' @examples
#' \donttest{
plotBox <- function(mr, x, y, ...)
{
params <- .extractParam(mr, x, y)
- L = length(params)
+ 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))
args <- list(...)
#' Note that the drawing does not correspond to a function; it is just a
#' convenient way to visualize the estimated parameters.
#'
+#' @name plotCoefs
+#'
#' @param mr List of parameters matrices
#' @param params True value of the parameters matrix
#' @param ... Additional graphical parameters
#' 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.
#'
+#' @name generateSampleIO
+#'
#' @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
#' io$index[1] #number of the group of X[1,] and Y[1] (in 1...K)
#'
#' @export
-generateSampleIO = function(n, p, β, b, link)
+generateSampleIO <- function(n, p, β, b, link)
{
# Check arguments
- tryCatch({n = as.integer(n)}, error=function(e) stop("Cannot convert n to integer"))
+ 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)
#' normalize(x) #column 1 is 1/sqrt(5) (1 2),
#' #and column 2 is 1/sqrt(10) (-1, 3)
#' @export
-normalize = function(x)
+normalize <- function(x)
{
- x = as.matrix(x)
- norm2 = sqrt( colSums(x^2) )
+ x <- as.matrix(x)
+ norm2 <- sqrt( colSums(x^2) )
sweep(x, 2, norm2, '/')
}
#
# @return Matrix of size dxd
#
-.T_I_I_w = function(Te, w)
+.T_I_I_w <- function(Te, w)
{
- d = length(w)
- Ma = matrix(0,nrow=d,ncol=d)
+ d <- length(w)
+ Ma <- matrix(0,nrow=d,ncol=d)
for (j in 1:d)
- Ma = Ma + w[j] * Te[,,j]
+ Ma <- Ma + w[j] * Te[,,j]
Ma
}
#
# @return Matrix of size dxd
#
-.Moments_M2 = function(X, Y)
+.Moments_M2 <- function(X, Y)
{
- n = nrow(X)
- d = ncol(X)
- M2 = matrix(0,nrow=d,ncol=d)
+ 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)
}
#
# @return Array of size dxdxd
#
-.Moments_M3 = function(X, Y)
+.Moments_M3 <- function(X, Y)
{
- n = nrow(X)
- d = ncol(X)
- M3 = array(0,dim=c(d,d,d))
+ 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) )
}
#
# @return A permutation minimizing cost
#
-.hungarianAlgorithm = function(distances)
+.hungarianAlgorithm <- function(distances)
{
- n = nrow(distances)
+ n <- nrow(distances)
.C("hungarianAlgorithm", distances=as.double(distances), pn=as.integer(n),
assignment=integer(n), PACKAGE="morpheus")$assignment
}
#' # a[[i]] is expected to contain m1 for all i
#'
#' @export
-alignMatrices = function(Ms, ref, ls_mode=c("exact","approx1","approx2"))
+alignMatrices <- function(Ms, ref, ls_mode=c("exact","approx1","approx2"))
{
if (!is.matrix(ref) || any(is.na(ref)))
stop("ref: matrix, no NAs")
L <- length(Ms)
for (i in 1:L)
{
- m = Ms[[i]] #shorthand
+ m <- Ms[[i]] #shorthand
if (ls_mode == "exact")
{
// Normalize W: x 1/n
for (int j=0; j<dim; j++)
{
- for (int k=0; k<=j; k++)
+ for (int k=j; k<dim; k++)
W[mi(j,k,dim,dim)] /= n;
}
- // Symmetrize W: W[j,k] = W[k,j] for k > j
+ // Symmetrize W: W[k,j] = W[j,k] for k > j
for (int j=0; j<dim; j++)
{
for (int k=j+1; k<dim; k++)
- W[mi(j,k,dim,dim)] = W[mi(k,j,dim,dim)];
+ W[mi(k,j,dim,dim)] = W[mi(j,k,dim,dim)];
}
free(g);
}
library(testthat)
-
-# Locally:
-#library(devtools)
-#load_all("../")
-
-# With R CMD check:
library(morpheus)
test_check("morpheus")
-context("Moments")
-
test_that("both versions of Moments_Mi agree on various inputs",
{
for (n in c(20,200))
-context("alignMatrices")
-
# Helper to generate a random series of matrices to align
.generateMatrices = function(d, K, N, noise)
{
# 3] Check alignment
for (j in 2:N)
- expect_that( norm(aligned[[j]] - Ms[[1]]), is_less_than(max_error) )
+ expect_lt( norm(aligned[[j]] - Ms[[1]]), max_error )
}
})
-context("computeMu")
-
test_that("on input of sufficient size, β/||β|| is estimated accurately enough",
{
n <- 100000
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) )
+ # NOTE: 0.5 is loose threshold, but values around 0.3 are expected...
+ expect_lt( diff_norm, 0.5 )
}
}
})
-context("hungarianAlgorithm")
-
test_that("HungarianAlgorithm provides the correct answer on various inputs",
{
for (n in c(2,3,10,50))
-context("jointDiag::ajd")
-
#auxiliary to test diagonality
.computeMuCheckDiag = function(X, Y, K, jd_method, β_ref)
{
for (i in 1:K)
{
shouldBeDiag <- invβ %*% M2_t[,,i] %*% t(invβ)
- expect_that(
+ expect_lt(
mean( abs(shouldBeDiag[upper.tri(shouldBeDiag) | lower.tri(shouldBeDiag)]) ),
- is_less_than(max_error) )
+ max_error)
}
}
-context("OptimParams")
-
naive_f <- function(link, M1,M2,M3, p,β,b)
{
d <- length(M1)
res
}
-# TODO: understand why it fails and reactivate this test
-#test_that("naive computation provides the same result as vectorized computations",
-#{
-# h <- 1e-7 #for finite-difference tests
-# tol <- 1e-3 #large tolerance, necessary in some cases... (generally 1e-6 is OK)
-# n <- 10
-# 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^2, -1, 1), ncol=d)
-# M3 <- array(runif(d^3, -1, 1), dim=c(d,d,d))
-#
-# for (link in c("logit","probit"))
-# {
-# # X and Y are unused here (W not re-computed)
-# op <- optimParams(X=matrix(runif(n*d),ncol=d), Y=rbinom(n,1,.5),
-# K, link, M=list(M1,M2,M3))
-# op$W <- diag(d + d^2 + d^3)
-#
-# 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)
-# }
-# }
-# }
-#})
+# TODO: understand why delta is so large (should be 10^-6 10^-7 ...)
+test_that("naive computation provides the same result as vectorized computations",
+{
+ h <- 1e-7 #for finite-difference tests
+ n <- 10
+ 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^2, -1, 1), ncol=d)
+ M3 <- array(runif(d^3, -1, 1), dim=c(d,d,d))
+
+ for (link in c("logit","probit"))
+ {
+ # X and Y are unused here (W not re-computed)
+ op <- optimParams(X=matrix(runif(n*d),ncol=d), Y=rbinom(n,1,.5),
+ K, link, M=list(M1,M2,M3))
+ op$W <- diag(d + d^2 + d^3)
+
+ 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 (TODO: 1 is way too high)
+ expect_equal( op$f(x)[1], naive_f(link,M1,M2,M3, p,β,b), tolerance=1 )
+
+ # 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)[1], tolerance=0.5 )
+ }
+ }
+ }
+})
test_that("W computed in C and in R are the same",
{
tol <- 1e-8
- n <- 500
- for (dK in list( c(2,2), c(5,3)))
+ n <- 10
+ for (dK in list( c(2,2))) #, c(5,3)))
{
d <- dK[1]
K <- dK[2]
W <- matrix(0, nrow=dd, ncol=dd)
Omega2 <- matrix( .C("Compute_Omega",
X=as.double(X), Y=as.integer(Y), M=as.double(M),
- pn=as.integer(n), pd=as.integer(d),
+ pnc=as.integer(1), pn=as.integer(n), pd=as.integer(d),
W=as.double(W), PACKAGE="morpheus")$W, nrow=dd, ncol=dd )
rg <- range(Omega1 - Omega2)
- expect_equal(rg[1], rg[2], tol)
+ expect_equal(rg[1], rg[2], tolerance=tol)
}
})
#!/bin/bash
N=100
-n=5e3
nc=3
#nstart=5
for n in "5000" "10000" "100000"; do
+ echo $n
for d in 2 5 10; do
- for link in "logit" "probit"; do
+ echo " "$d
+ for link in "logit"; do #"probit"; do
#R --slave --args N=$N n=$n nc=$nc d=$d link=$link nstart=$nstart <multistart.R >out_${n}_${link}_${d}_${nstart} 2>&1
R --slave --args N=$N n=$n nc=$nc d=$d link=$link <accuracy.R >out_${n}_${link}_${d} 2>&1
done
projects / morpheus.git / commitdiff