X-Git-Url: https://git.auder.net/?p=valse.git;a=blobdiff_plain;f=pkg%2FR%2FEMGLLF.R;h=2aeea537b9d6656c04032c2697d407f6c4d1d823;hp=f944f98e38ca48fcac75c73cd3ac2074551d06e9;hb=17b9fa5f6004bb55e915d8916f1c9a1f128a65ce;hpb=ca277ac5ab51fef149014eb5e4610403fdb3227b diff --git a/pkg/R/EMGLLF.R b/pkg/R/EMGLLF.R index f944f98..2aeea53 100644 --- a/pkg/R/EMGLLF.R +++ b/pkg/R/EMGLLF.R @@ -19,7 +19,8 @@ #' rho : parametre de variance renormalisé, calculé par l'EM #' pi : parametre des proportions renormalisé, calculé par l'EM #' LLF : log vraisemblance associée à cet échantillon, pour les valeurs estimées des paramètres -#' S : ... affec : ... +#' S : ... +#' affec : ... #' #' @export EMGLLF <- function(phiInit, rhoInit, piInit, gamInit, mini, maxi, gamma, lambda, @@ -41,21 +42,27 @@ EMGLLF <- function(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, PACKAGE = "valse") + list(phi = phi, rho = rho, pi = pi, llh = llh, S = S, affec=affec) } # R version - slow but easy to read .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 +162,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) { @@ -185,5 +192,6 @@ EMGLLF <- function(phiInit, rhoInit, piInit, gamInit, mini, maxi, gamma, lambda, break } - list(phi = phi, rho = rho, pi = pi, llh = llh, S = S) + affec = apply(gam, 1, which.max) + list(phi = phi, rho = rho, pi = pi, llh = llh, S = S, affec=affec) }