X-Git-Url: https://git.auder.net/?p=valse.git;a=blobdiff_plain;f=pkg%2FR%2FinitSmallEM.R;h=bd5ce1779303f3d17bf727b30048d5768b7a12a6;hp=541d7e1cb922b590d1893eea5bba3472eeb1e4ed;hb=4d9db27f0d1749e5577038dedbc5f4d0826f2772;hpb=f33f35efc9a01f93bb61959522d90ee6a76b892e diff --git a/pkg/R/initSmallEM.R b/pkg/R/initSmallEM.R index 541d7e1..bd5ce17 100644 --- a/pkg/R/initSmallEM.R +++ b/pkg/R/initSmallEM.R @@ -8,26 +8,26 @@ #' @export #' @importFrom methods new #' @importFrom stats cutree dist hclust runif -initSmallEM = function(k,X,Y) +initSmallEM = function(k,X,Y, fast=TRUE) { n = nrow(Y) m = ncol(Y) p = ncol(X) - - Zinit1 = array(0, dim=c(n,20)) - betaInit1 = array(0, dim=c(p,m,k,20)) - sigmaInit1 = array(0, dim = c(m,m,k,20)) - phiInit1 = array(0, dim = c(p,m,k,20)) - rhoInit1 = array(0, dim = c(m,m,k,20)) + nIte = 20 + Zinit1 = array(0, dim=c(n,nIte)) + betaInit1 = array(0, dim=c(p,m,k,nIte)) + sigmaInit1 = array(0, dim = c(m,m,k,nIte)) + phiInit1 = array(0, dim = c(p,m,k,nIte)) + rhoInit1 = array(0, dim = c(m,m,k,nIte)) Gam = matrix(0, n, k) - piInit1 = matrix(0,20,k) - gamInit1 = array(0, dim=c(n,k,20)) + 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:20) + #require(MASS) #Moore-Penrose generalized inverse of matrix + for(repet in 1:nIte) { - distance_clus = dist(X) + distance_clus = dist(cbind(X,Y)) tree_hier = hclust(distance_clus) Zinit1[,repet] = cutree(tree_hier, k) @@ -36,10 +36,10 @@ initSmallEM = function(k,X,Y) Z = Zinit1[,repet] Z_indice = seq_len(n)[Z == r] #renvoit les indices où Z==r if (length(Z_indice) == 1) { - betaInit1[,,r,repet] = 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] = 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) @@ -62,14 +62,12 @@ initSmallEM = function(k,X,Y) miniInit = 10 maxiInit = 11 - #new_EMG = .Call("EMGLLF_core",phiInit1[,,,repet],rhoInit1[,,,repet],piInit1[repet,], -# gamInit1[,,repet],miniInit,maxiInit,1,0,X,Y,1e-4) - new_EMG = EMGLLF(phiInit1[,,,repet],rhoInit1[,,,repet],piInit1[repet,],gamInit1[,,repet],miniInit,maxiInit,1,0,X,Y,1e-4) - LLFEessai = new_EMG$LLF + init_EMG = EMGLLF(phiInit1[,,,repet], rhoInit1[,,,repet], piInit1[repet,], + gamInit1[,,repet], miniInit, maxiInit, gamma=1, lambda=0, X, Y, eps=1e-4, fast) + LLFEessai = init_EMG$LLF LLFinit1[repet] = LLFEessai[length(LLFEessai)] } - - b = which.max(LLFinit1) + b = which.min(LLFinit1) phiInit = phiInit1[,,,b] rhoInit = rhoInit1[,,,b] piInit = piInit1[b,]