X-Git-Url: https://git.auder.net/?p=valse.git;a=blobdiff_plain;f=pkg%2FR%2FinitSmallEM.R;h=44b4b06afa909893bf44a91f8c699af4cbb9aa99;hp=fafa2c464ce5a16a87760e9d67374e9a23278cf7;hb=ea5860f1b4fc91f06e371a0b26915198474a849d;hpb=ffdf94474d96cdd3e9d304ce809df7e62aa957ed diff --git a/pkg/R/initSmallEM.R b/pkg/R/initSmallEM.R index fafa2c4..44b4b06 100644 --- a/pkg/R/initSmallEM.R +++ b/pkg/R/initSmallEM.R @@ -8,11 +8,11 @@ #' @export #' @importFrom methods new #' @importFrom stats cutree dist hclust runif -initSmallEM <- function(k, X, Y, fast = TRUE) +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)) @@ -23,24 +23,22 @@ initSmallEM <- function(k, X, Y, fast = TRUE) piInit1 <- matrix(0, nIte, k) gamInit1 <- array(0, dim = c(n, k, nIte)) LLFinit1 <- list() - + # require(MASS) #Moore-Penrose generalized inverse of matrix for (repet in 1:nIte) { distance_clus <- dist(cbind(X, Y)) tree_hier <- hclust(distance_clus) Zinit1[, repet] <- cutree(tree_hier, k) - + for (r in 1:k) { Z <- Zinit1[, repet] Z_indice <- seq_len(n)[Z == r] #renvoit les indices où Z==r - if (length(Z_indice) == 1) - { + if (length(Z_indice) == 1) { betaInit1[, , r, repet] <- MASS::ginv(crossprod(t(X[Z_indice, ]))) %*% crossprod(t(X[Z_indice, ]), Y[Z_indice, ]) - } else - { + } else { betaInit1[, , r, repet] <- MASS::ginv(crossprod(X[Z_indice, ])) %*% crossprod(X[Z_indice, ], Y[Z_indice, ]) } @@ -49,34 +47,33 @@ initSmallEM <- function(k, X, Y, fast = TRUE) rhoInit1[, , r, repet] <- solve(sigmaInit1[, , r, repet]) piInit1[repet, r] <- mean(Z == r) } - + for (i in 1:n) { for (r in 1:k) { - 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) + dotProduct <- tcrossprod(Y[i, ] %*% rhoInit1[, , r, repet] + - X[i, ] %*% phiInit1[, , r, repet]) + Gam[i, r] <- piInit1[repet, r] * + gdet(rhoInit1[, , r, repet]) * exp(-0.5 * dotProduct) } sumGamI <- sum(Gam[i, ]) gamInit1[i, , repet] <- Gam[i, ]/sumGamI } - + miniInit <- 10 maxiInit <- 11 - - init_EMG <- EMGLLF(phiInit1[, , , repet], rhoInit1[, , , repet], piInit1[repet, - ], gamInit1[, , repet], miniInit, maxiInit, gamma = 1, lambda = 0, X, - Y, eps = 1e-04, fast) - LLFEessai <- init_EMG$LLF - LLFinit1[repet] <- LLFEessai[length(LLFEessai)] + + init_EMG <- EMGLLF(phiInit1[, , , repet], rhoInit1[, , , repet], piInit1[repet, ], + gamInit1[, , repet], miniInit, maxiInit, gamma = 1, lambda = 0, X, Y, + eps = 1e-04, fast) + LLFinit1[[repet]] <- init_EMG$llh } b <- which.min(LLFinit1) phiInit <- phiInit1[, , , b] rhoInit <- rhoInit1[, , , b] piInit <- piInit1[b, ] gamInit <- gamInit1[, , b] - + return(list(phiInit = phiInit, rhoInit = rhoInit, piInit = piInit, gamInit = gamInit)) }