X-Git-Url: https://git.auder.net/?p=valse.git;a=blobdiff_plain;f=pkg%2FR%2FEMGrank.R;h=db409481bcfffe6a94000d6982b18ce1f2ccedb1;hp=436b28982eb093c3c4eda822399ff50b850ffa29;hb=1b698c1619dbcf5b3a0608dc894d249945d2bce3;hpb=f7e157cdbcf2d60224c2d6773da9c698174e9aee diff --git a/pkg/R/EMGrank.R b/pkg/R/EMGrank.R index 436b289..db40948 100644 --- a/pkg/R/EMGrank.R +++ b/pkg/R/EMGrank.R @@ -1,4 +1,4 @@ -#' EMGrank +#' EMGrank #' #' Description de EMGrank #' @@ -23,7 +23,7 @@ EMGrank <- function(Pi, Rho, mini, maxi, X, Y, tau, rank, fast = TRUE) # Function in R return(.EMGrank_R(Pi, Rho, mini, maxi, X, Y, tau, rank)) } - + # Function in C n <- nrow(X) #nombre d'echantillons p <- ncol(X) #nombre de covariables @@ -50,12 +50,12 @@ matricize <- function(X) p <- dim(X)[2] m <- dim(Rho)[2] k <- dim(Rho)[3] - + # init outputs phi <- array(0, dim = c(p, m, k)) Z <- rep(1, n) LLF <- 0 - + # local variables Phi <- array(0, dim = c(p, m, k)) deltaPhi <- c() @@ -69,12 +69,12 @@ matricize <- function(X) # M step: update for Beta ( and then phi) for (r in 1:k) { - Z_indice <- seq_len(n)[Z == r] #indices where Z == r + Z_indice <- seq_len(n)[Z == r] #indices where Z == r if (length(Z_indice) == 0) next # U,S,V = SVD of (t(Xr)Xr)^{-1} * t(Xr) * Yr - s <- svd(MASS::ginv(crossprod(matricize(X[Z_indice, ]))) %*% crossprod(matricize(X[Z_indice, - ]), matricize(Y[Z_indice, ]))) + s <- svd(MASS::ginv(crossprod(matricize(X[Z_indice, ]))) + %*% crossprod(matricize(X[Z_indice, ]), matricize(Y[Z_indice, ]))) S <- s$d # Set m-rank(r) singular values to zero, and recompose best rank(r) approximation # of the initial product @@ -82,7 +82,7 @@ matricize <- function(X) S[(rank[r] + 1):length(S)] <- 0 phi[, , r] <- s$u %*% diag(S) %*% t(s$v) %*% Rho[, , r] } - + # Step E and computation of the loglikelihood sumLogLLF2 <- 0 for (i in seq_len(n)) @@ -91,8 +91,7 @@ matricize <- function(X) maxLogGamIR <- -Inf for (r in seq_len(k)) { - dotProduct <- tcrossprod(Y[i, ] %*% Rho[, , r] - X[i, ] %*% phi[, - , r]) + dotProduct <- tcrossprod(Y[i, ] %*% Rho[, , r] - X[i, ] %*% phi[, , r]) logGamIR <- log(Pi[r]) + log(det(Rho[, , r])) - 0.5 * dotProduct # Z[i] = index of max (gam[i,]) if (logGamIR > maxLogGamIR) @@ -104,15 +103,15 @@ matricize <- function(X) } sumLogLLF2 <- sumLogLLF2 + log(sumLLF1) } - + LLF <- -1/n * sumLogLLF2 - + # update distance parameter to check algorithm convergence (delta(phi, Phi)) - deltaPhi <- c(deltaPhi, max((abs(phi - Phi))/(1 + abs(phi)))) #TODO: explain? + deltaPhi <- c(deltaPhi, max((abs(phi - Phi))/(1 + abs(phi)))) #TODO: explain? if (length(deltaPhi) > deltaPhiBufferSize) deltaPhi <- deltaPhi[2:length(deltaPhi)] sumDeltaPhi <- sum(abs(deltaPhi)) - + # update other local variables Phi <- phi ite <- ite + 1