X-Git-Url: https://git.auder.net/?p=valse.git;a=blobdiff_plain;f=pkg%2FR%2FinitSmallEM.R;h=39453227c5847f0b4341de18287ae5d773492e3b;hp=056d7e70e89fb08e7072ac63820c9c5f6365dc8a;hb=6af1d4897dbab92a7be05068e0e15823378965d9;hpb=a3cbbaea1cc3c107e5ca62ed1ffe7b9499de0a91 diff --git a/pkg/R/initSmallEM.R b/pkg/R/initSmallEM.R index 056d7e7..3945322 100644 --- a/pkg/R/initSmallEM.R +++ b/pkg/R/initSmallEM.R @@ -1,18 +1,24 @@ -#' initialization of the EM algorithm +#' initSmallEM +#' +#' initialization of the EM algorithm #' #' @param k number of components #' @param X matrix of covariates (of size n*p) #' @param Y matrix of responses (of size n*m) +#' @param fast boolean to enable or not the C function call +#' +#' @return a list with phiInit (the regression parameter reparametrized), +#' rhoInit (the covariance parameter reparametrized), piInit (the proportion parameter is the +#' mixture model), gamInit (the conditional expectation) #' -#' @return a list with phiInit, rhoInit, piInit, gamInit -#' @export -#' @importFrom methods new #' @importFrom stats cutree dist hclust runif +#' +#' @export 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)) @@ -34,12 +40,12 @@ initSmallEM <- function(k, X, Y, fast) for (r in 1:k) { Z <- Zinit1[, repet] - Z_indice <- seq_len(n)[Z == r] #renvoit les indices où Z==r + Z_indice <- seq_len(n)[Z == r] #renvoit les indices ou Z==r if (length(Z_indice) == 1) { - betaInit1[, , r, repet] <- MASS::ginv(crossprod(t(X[Z_indice, ]))) %*% + betaInit1[, , r, repet] <- MASS::ginv(crossprod(t(X[Z_indice, ]))) %*% crossprod(t(X[Z_indice, ]), Y[Z_indice, ]) } else { - betaInit1[, , r, repet] <- MASS::ginv(crossprod(X[Z_indice, ])) %*% + betaInit1[, , r, repet] <- MASS::ginv(crossprod(X[Z_indice, ])) %*% crossprod(X[Z_indice, ], Y[Z_indice, ]) } sigmaInit1[, , r, repet] <- diag(m) @@ -54,10 +60,11 @@ initSmallEM <- function(k, X, Y, fast) { dotProduct <- tcrossprod(Y[i, ] %*% rhoInit1[, , r, repet] - X[i, ] %*% phiInit1[, , r, repet]) - Gam[i, r] <- piInit1[repet, r] * + Gam[i, r] <- piInit1[repet, r] * det(rhoInit1[, , r, repet]) * exp(-0.5 * dotProduct) } sumGamI <- sum(Gam[i, ]) + # TODO: next line is a division by zero if dotProduct is big gamInit1[i, , repet] <- Gam[i, ]/sumGamI } @@ -75,5 +82,5 @@ initSmallEM <- function(k, X, Y, fast) piInit <- piInit1[b, ] gamInit <- gamInit1[, , b] - return(list(phiInit = phiInit, rhoInit = rhoInit, piInit = piInit, gamInit = gamInit)) + list(phiInit = phiInit, rhoInit = rhoInit, piInit = piInit, gamInit = gamInit) }