From: Benjamin Auder Date: Fri, 3 Nov 2017 09:54:04 +0000 (+0100) Subject: progress in debug: fix LLF/llh mismatch, and length + add adapter test in test/ X-Git-Url: https://git.auder.net/?p=valse.git;a=commitdiff_plain;h=6279ba8656582370e7242ff9e77d22c23fe8ca04 progress in debug: fix LLF/llh mismatch, and length + add adapter test in test/ --- diff --git a/pkg/R/EMGLLF.R b/pkg/R/EMGLLF.R index 08ff203..57638f9 100644 --- a/pkg/R/EMGLLF.R +++ b/pkg/R/EMGLLF.R @@ -40,7 +40,7 @@ EMGLLF <- function(phiInit, rhoInit, piInit, gamInit, mini, maxi, gamma, lambda, k <- length(piInit) #nombre de composantes dans le mélange .Call("EMGLLF", phiInit, rhoInit, piInit, gamInit, mini, maxi, gamma, lambda, X, Y, eps, phi = double(p * m * k), rho = double(m * m * k), pi = double(k), - LLF = double(maxi), S = double(p * m * k), affec = integer(n), n, p, m, k, + llh = double(1), S = double(p * m * k), affec = integer(n), n, p, m, k, PACKAGE = "valse") } diff --git a/pkg/R/EMGrank.R b/pkg/R/EMGrank.R index b85a0fa..f2bf58e 100644 --- a/pkg/R/EMGrank.R +++ b/pkg/R/EMGrank.R @@ -8,7 +8,7 @@ #' @param maxi Nombre maximal d'itérations dans l'algorithme EM #' @param X Régresseurs #' @param Y Réponse -#' @param tau Seuil pour accepter la convergence +#' @param eps Seuil pour accepter la convergence #' @param rank Vecteur des rangs possibles #' #' @return A list ... @@ -16,12 +16,12 @@ #' LLF : log vraisemblance associé à cet échantillon, pour les valeurs estimées des paramètres #' #' @export -EMGrank <- function(Pi, Rho, mini, maxi, X, Y, tau, rank, fast = TRUE) +EMGrank <- function(Pi, Rho, mini, maxi, X, Y, eps, rank, fast = TRUE) { if (!fast) { # Function in R - return(.EMGrank_R(Pi, Rho, mini, maxi, X, Y, tau, rank)) + return(.EMGrank_R(Pi, Rho, mini, maxi, X, Y, eps, rank)) } # Function in C @@ -29,7 +29,7 @@ EMGrank <- function(Pi, Rho, mini, maxi, X, Y, tau, rank, fast = TRUE) p <- ncol(X) #nombre de covariables m <- ncol(Y) #taille de Y (multivarié) k <- length(Pi) #nombre de composantes dans le mélange - .Call("EMGrank", Pi, Rho, mini, maxi, X, Y, tau, rank, phi = double(p * m * k), + .Call("EMGrank", Pi, Rho, mini, maxi, X, Y, eps, rank, phi = double(p * m * k), LLF = double(1), n, p, m, k, PACKAGE = "valse") } @@ -43,7 +43,7 @@ matricize <- function(X) } # R version - slow but easy to read -.EMGrank_R <- function(Pi, Rho, mini, maxi, X, Y, tau, rank) +.EMGrank_R <- function(Pi, Rho, mini, maxi, X, Y, eps, rank) { # matrix dimensions n <- nrow(X) @@ -64,7 +64,7 @@ matricize <- function(X) # main loop ite <- 1 - while (ite <= mini || (ite <= maxi && sumDeltaPhi > tau)) + while (ite <= mini || (ite <= maxi && sumDeltaPhi > eps)) { # M step: update for Beta ( and then phi) for (r in 1:k) diff --git a/pkg/R/computeGridLambda.R b/pkg/R/computeGridLambda.R index 8449d10..ac0788a 100644 --- a/pkg/R/computeGridLambda.R +++ b/pkg/R/computeGridLambda.R @@ -11,13 +11,13 @@ #' @param gamma power of weights in the penalty #' @param mini minimum number of iterations in EM algorithm #' @param maxi maximum number of iterations in EM algorithm -#' @param tau threshold to stop EM algorithm +#' @param eps threshold to stop EM algorithm #' #' @return the grid of regularization parameters #' #' @export computeGridLambda <- function(phiInit, rhoInit, piInit, gamInit, X, Y, gamma, mini, - maxi, tau, fast) + maxi, eps, fast) { n <- nrow(X) p <- ncol(X) @@ -25,7 +25,8 @@ computeGridLambda <- function(phiInit, rhoInit, piInit, gamInit, X, Y, gamma, mi k <- length(piInit) list_EMG <- EMGLLF(phiInit, rhoInit, piInit, gamInit, mini, maxi, gamma, lambda = 0, - X, Y, tau, fast) + X, Y, eps, fast) + grid <- array(0, dim = c(p, m, k)) for (j in 1:p) { diff --git a/pkg/R/main.R b/pkg/R/main.R index d29fe69..387d553 100644 --- a/pkg/R/main.R +++ b/pkg/R/main.R @@ -123,8 +123,7 @@ valse <- function(X, Y, procedure = "LassoMLE", selecMod = "DDSE", gamma = 1, mi complexity = sumPen, contrast = -LLH) })) tableauRecap <- tableauRecap[which(tableauRecap[, 4] != Inf), ] - - + if (verbose == TRUE) print(tableauRecap) modSel <- capushe::capushe(tableauRecap, n) @@ -141,11 +140,11 @@ valse <- function(X, Y, procedure = "LassoMLE", selecMod = "DDSE", gamma = 1, mi { modSel@AIC_capushe$model } - + listMod <- as.integer(unlist(strsplit(as.character(indModSel), "[.]"))) modelSel <- models_list[[listMod[1]]][[listMod[2]]] modelSel$tableau <- tableauRecap - + if (plot) print(plot_valse(X, Y, modelSel, n)) diff --git a/pkg/src/adapters/a.EMGLLF.c b/pkg/src/adapters/a.EMGLLF.c index e93b7ec..9b004c2 100644 --- a/pkg/src/adapters/a.EMGLLF.c +++ b/pkg/src/adapters/a.EMGLLF.c @@ -46,7 +46,7 @@ SEXP EMGLLF( // OUTPUTS // ///////////// - SEXP phi, rho, pi, LLF, S, affec, dimPhiS, dimRho; + SEXP phi, rho, pi, llh, S, affec, dimPhiS, dimRho; PROTECT(dimPhiS = allocVector(INTSXP, 3)); int* pDimPhiS = INTEGER(dimPhiS); pDimPhiS[0] = p; pDimPhiS[1] = m; pDimPhiS[2] = k; @@ -56,10 +56,10 @@ SEXP EMGLLF( PROTECT(phi = allocArray(REALSXP, dimPhiS)); PROTECT(rho = allocArray(REALSXP, dimRho)); PROTECT(pi = allocVector(REALSXP, k)); - PROTECT(LLF = allocVector(REALSXP, maxi)); + PROTECT(llh = allocVector(REALSXP, 1)); PROTECT(S = allocArray(REALSXP, dimPhiS)); PROTECT(affec = allocVector(INTSXP, n)); - double *pPhi=REAL(phi), *pRho=REAL(rho), *pPi=REAL(pi), *pLLF=REAL(LLF), *pS=REAL(S); + double *pPhi=REAL(phi), *pRho=REAL(rho), *pPi=REAL(pi), *pLlh=REAL(llh), *pS=REAL(S); int *pAffec=INTEGER(affec); //////////////////// @@ -67,14 +67,14 @@ SEXP EMGLLF( //////////////////// EMGLLF_core(phiInit,rhoInit,piInit,gamInit,mini,maxi,gamma,lambda,X,Y,eps, - pPhi,pRho,pPi,pLLF,pS,pAffec, + pPhi,pRho,pPi,pLlh,pS,pAffec, n,p,m,k); // Build list from OUT params and return it SEXP listParams, listNames; int nouts = 6; PROTECT(listParams = allocVector(VECSXP, nouts)); - char* lnames[6] = {"phi", "rho", "pi", "LLF", "S", "affec"}; //lists labels + char* lnames[6] = {"phi", "rho", "pi", "llh", "S", "affec"}; //lists labels PROTECT(listNames = allocVector(STRSXP,nouts)); for (int i=0; i