From ea5860f1b4fc91f06e371a0b26915198474a849d Mon Sep 17 00:00:00 2001 From: Benjamin Auder Date: Fri, 21 Apr 2017 15:33:07 +0200 Subject: [PATCH] fix for m==1 --- pkg/DESCRIPTION | 1 + pkg/R/EMGLLF.R | 21 +++++++++++++-------- pkg/R/EMGrank.R | 10 +++++----- pkg/R/computeGridLambda.R | 6 +++--- pkg/R/constructionModelesLassoMLE.R | 14 +++++++------- pkg/R/constructionModelesLassoRank.R | 6 +++--- pkg/R/initSmallEM.R | 6 +++--- pkg/R/main.R | 8 ++++---- pkg/R/selectVariables.R | 13 ++++++++++--- pkg/R/util.R | 7 +++++++ 10 files changed, 56 insertions(+), 36 deletions(-) create mode 100644 pkg/R/util.R diff --git a/pkg/DESCRIPTION b/pkg/DESCRIPTION index e18fddb..3b33e25 100644 --- a/pkg/DESCRIPTION +++ b/pkg/DESCRIPTION @@ -39,3 +39,4 @@ Collate: 'EMGLLF.R' 'generateXY.R' 'A_NAMESPACE.R' + 'util.R' diff --git a/pkg/R/EMGLLF.R b/pkg/R/EMGLLF.R index 6ee7ba7..03f0a75 100644 --- a/pkg/R/EMGLLF.R +++ b/pkg/R/EMGLLF.R @@ -47,15 +47,20 @@ EMGLLF <- function(phiInit, rhoInit, piInit, gamInit, mini, maxi, gamma, lambda, .EMGLLF_R <- function(phiInit, rhoInit, piInit, gamInit, mini, maxi, gamma, lambda, X, Y, eps) { - # Matrix dimensions: NOTE: phiInit *must* be an array (even if p==1) - n <- dim(Y)[1] - p <- dim(phiInit)[1] - m <- dim(phiInit)[2] - k <- dim(phiInit)[3] + # Matrix dimensions + n <- nrow(X) + p <- ncol(X) + m <- ncol(Y) + k <- length(piInit) + + # Adjustments required when p==1 or m==1 (var.sel. or output dim 1) + if (p==1 || m==1) + phiInit <- array(phiInit, dim=c(p,m,k)) + if (m==1) + rhoInit <- array(rhoInit, dim=c(m,m,k)) # Outputs - phi <- array(NA, dim = c(p, m, k)) - phi[1:p, , ] <- phiInit + phi <- phiInit rho <- rhoInit pi <- piInit llh <- -Inf @@ -155,7 +160,7 @@ EMGLLF <- function(phiInit, rhoInit, piInit, gamInit, mini, maxi, gamma, lambda, ## E step # Precompute det(rho[,,r]) for r in 1...k - detRho <- sapply(1:k, function(r) det(rho[, , r])) + detRho <- sapply(1:k, function(r) gdet(rho[, , r])) sumLogLLH <- 0 for (i in 1:n) { diff --git a/pkg/R/EMGrank.R b/pkg/R/EMGrank.R index db0b8f2..b85a0fa 100644 --- a/pkg/R/EMGrank.R +++ b/pkg/R/EMGrank.R @@ -46,10 +46,10 @@ matricize <- function(X) .EMGrank_R <- function(Pi, Rho, mini, maxi, X, Y, tau, rank) { # matrix dimensions - n <- dim(X)[1] - p <- dim(X)[2] - m <- dim(Rho)[2] - k <- dim(Rho)[3] + n <- nrow(X) + p <- ncol(X) + m <- ncol(Y) + k <- length(Pi) # init outputs phi <- array(0, dim = c(p, m, k)) @@ -92,7 +92,7 @@ matricize <- function(X) for (r in seq_len(k)) { dotProduct <- tcrossprod(Y[i, ] %*% Rho[, , r] - X[i, ] %*% phi[, , r]) - logGamIR <- log(Pi[r]) + log(det(Rho[, , r])) - 0.5 * dotProduct + logGamIR <- log(Pi[r]) + log(gdet(Rho[, , r])) - 0.5 * dotProduct # Z[i] = index of max (gam[i,]) if (logGamIR > maxLogGamIR) { diff --git a/pkg/R/computeGridLambda.R b/pkg/R/computeGridLambda.R index c2e9c8c..cf762ec 100644 --- a/pkg/R/computeGridLambda.R +++ b/pkg/R/computeGridLambda.R @@ -20,9 +20,9 @@ computeGridLambda <- function(phiInit, rhoInit, piInit, gamInit, X, Y, gamma, mi maxi, tau, fast) { n <- nrow(X) - p <- dim(phiInit)[1] - m <- dim(phiInit)[2] - k <- dim(phiInit)[3] + p <- ncol(X) + m <- ncol(Y) + k <- length(piInit) list_EMG <- EMGLLF(phiInit, rhoInit, piInit, gamInit, mini, maxi, gamma, lambda = 0, X, Y, tau, fast) diff --git a/pkg/R/constructionModelesLassoMLE.R b/pkg/R/constructionModelesLassoMLE.R index 90d0a2a..1275ca3 100644 --- a/pkg/R/constructionModelesLassoMLE.R +++ b/pkg/R/constructionModelesLassoMLE.R @@ -40,10 +40,10 @@ constructionModelesLassoMLE <- function(phiInit, rhoInit, piInit, gamInit, mini, if (verbose) print(paste("Computations for lambda=", lambda)) - n <- dim(X)[1] - p <- dim(phiInit)[1] - m <- dim(phiInit)[2] - k <- dim(phiInit)[3] + n <- nrow(X) + p <- ncol(X) + m <- ncol(Y) + k <- length(piInit) sel.lambda <- S[[lambda]]$selected # col.sel = which(colSums(sel.lambda)!=0) #if boolean matrix col.sel <- which(sapply(sel.lambda, length) > 0) #if list of selected vars @@ -51,8 +51,8 @@ constructionModelesLassoMLE <- function(phiInit, rhoInit, piInit, gamInit, mini, return(NULL) # lambda == 0 because we compute the EMV: no penalization here - res <- EMGLLF(array(phiInit[col.sel, , ],dim=c(length(col.sel),m,k)), rhoInit, - piInit, gamInit, mini, maxi, gamma, 0, as.matrix(X[, col.sel]), Y, eps, fast) + res <- EMGLLF(array(phiInit,dim=c(p,m,k))[col.sel, , ], rhoInit, piInit, gamInit, + mini, maxi, gamma, 0, as.matrix(X[, col.sel]), Y, eps, fast) # Eval dimension from the result + selected phiLambda2 <- res$phi @@ -71,7 +71,7 @@ constructionModelesLassoMLE <- function(phiInit, rhoInit, piInit, gamInit, mini, { delta <- (Y %*% rhoLambda[, , r] - (X[, col.sel] %*% t(phiLambda[col.sel, , r]))) } else delta <- (Y %*% rhoLambda[, , r] - (X[, col.sel] %*% phiLambda[col.sel, , r])) - densite <- densite + piLambda[r] * det(rhoLambda[, , r])/(sqrt(2 * base::pi))^m * + densite <- densite + piLambda[r] * gdet(rhoLambda[, , r])/(sqrt(2 * base::pi))^m * exp(-diag(tcrossprod(delta))/2) } llhLambda <- c(sum(log(densite)), (dimension + m + 1) * k - 1) diff --git a/pkg/R/constructionModelesLassoRank.R b/pkg/R/constructionModelesLassoRank.R index 85685e9..dc88f67 100644 --- a/pkg/R/constructionModelesLassoRank.R +++ b/pkg/R/constructionModelesLassoRank.R @@ -21,9 +21,9 @@ constructionModelesLassoRank <- function(S, k, mini, maxi, X, Y, eps, rank.min, rank.max, ncores, fast, verbose) { - n <- dim(X)[1] - p <- dim(X)[2] - m <- dim(Y)[2] + n <- nrow(X) + p <- ncol(X) + m <- ncol(Y) L <- length(S) # Possible interesting ranks diff --git a/pkg/R/initSmallEM.R b/pkg/R/initSmallEM.R index 056d7e7..44b4b06 100644 --- a/pkg/R/initSmallEM.R +++ b/pkg/R/initSmallEM.R @@ -10,9 +10,9 @@ #' @importFrom stats cutree dist hclust runif initSmallEM <- function(k, X, Y, fast) { - n <- nrow(Y) - m <- ncol(Y) + n <- nrow(X) p <- ncol(X) + m <- ncol(Y) nIte <- 20 Zinit1 <- array(0, dim = c(n, nIte)) betaInit1 <- array(0, dim = c(p, m, k, nIte)) @@ -55,7 +55,7 @@ initSmallEM <- function(k, X, Y, fast) dotProduct <- tcrossprod(Y[i, ] %*% rhoInit1[, , r, repet] - X[i, ] %*% phiInit1[, , r, repet]) Gam[i, r] <- piInit1[repet, r] * - det(rhoInit1[, , r, repet]) * exp(-0.5 * dotProduct) + gdet(rhoInit1[, , r, repet]) * exp(-0.5 * dotProduct) } sumGamI <- sum(Gam[i, ]) gamInit1[i, , repet] <- Gam[i, ]/sumGamI diff --git a/pkg/R/main.R b/pkg/R/main.R index fecf519..e741d65 100644 --- a/pkg/R/main.R +++ b/pkg/R/main.R @@ -31,9 +31,9 @@ valse <- function(X, Y, procedure = "LassoMLE", selecMod = "DDSE", gamma = 1, mi ncores_inner = 1, thresh = 1e-08, size_coll_mod = 10, fast = TRUE, verbose = FALSE, plot = TRUE) { - p <- dim(X)[2] - m <- dim(Y)[2] - n <- dim(X)[1] + n <- nrow(X) + p <- ncol(X) + m <- ncol(Y) if (verbose) print("main loop: over all k and all lambda") @@ -138,7 +138,7 @@ valse <- function(X, Y, procedure = "LassoMLE", selecMod = "DDSE", gamma = 1, mi for (r in 1:length(modelSel$pi)) { sqNorm2 <- sum((Y[i, ] %*% modelSel$rho[, , r] - X[i, ] %*% modelSel$phi[, , r])^2) - Gam[i, r] <- modelSel$pi[r] * exp(-0.5 * sqNorm2) * det(modelSel$rho[, , r]) + Gam[i, r] <- modelSel$pi[r] * exp(-0.5 * sqNorm2) * gdet(modelSel$rho[, , r]) } } Gam <- Gam/rowSums(Gam) diff --git a/pkg/R/selectVariables.R b/pkg/R/selectVariables.R index bfe4042..d863a4b 100644 --- a/pkg/R/selectVariables.R +++ b/pkg/R/selectVariables.R @@ -37,15 +37,22 @@ selectVariables <- function(phiInit, rhoInit, piInit, gamInit, mini, maxi, gamma params <- EMGLLF(phiInit, rhoInit, piInit, gamInit, mini, maxi, gamma, lambda, X, Y, eps, fast) - p <- dim(phiInit)[1] - m <- dim(phiInit)[2] + p <- ncol(X) + m <- ncol(Y) # selectedVariables: list where element j contains vector of selected variables # in [1,m] selectedVariables <- lapply(1:p, function(j) { # from boolean matrix mxk of selected variables obtain the corresponding boolean # m-vector, and finally return the corresponding indices - seq_len(m)[apply(abs(params$phi[j, , ]) > thresh, 1, any)] + if (m>1) { + seq_len(m)[apply(abs(params$phi[j, , ]) > thresh, 1, any)] + } else { + if (any(params$phi[j, , ] > thresh)) + 1 + else + numeric(0) + } }) list(selected = selectedVariables, Rho = params$rho, Pi = params$pi) diff --git a/pkg/R/util.R b/pkg/R/util.R new file mode 100644 index 0000000..f8b01cc --- /dev/null +++ b/pkg/R/util.R @@ -0,0 +1,7 @@ +# ... +gdet <- function(M) +{ + if (is.matrix(M)) + return (det(M)) + return (M[1]) #numeric, double +} -- 2.44.0