X-Git-Url: https://git.auder.net/?p=valse.git;a=blobdiff_plain;f=pkg%2FR%2FinitSmallEM.R;h=056d7e70e89fb08e7072ac63820c9c5f6365dc8a;hp=541d7e1cb922b590d1893eea5bba3472eeb1e4ed;hb=a3cbbaea1cc3c107e5ca62ed1ffe7b9499de0a91;hpb=f87ff0f5116c0c1c59c5608e46563ff0f79e5d43

diff --git a/pkg/R/initSmallEM.R b/pkg/R/initSmallEM.R
index 541d7e1..056d7e7 100644
--- a/pkg/R/initSmallEM.R
+++ b/pkg/R/initSmallEM.R
@@ -1,4 +1,4 @@
-#' initialization of the EM algorithm
+#' initialization of the EM algorithm 
 #'
 #' @param k number of components
 #' @param X matrix of covariates (of size n*p)
@@ -8,72 +8,72 @@
 #' @export
 #' @importFrom methods new
 #' @importFrom stats cutree dist hclust runif
-initSmallEM = function(k,X,Y)
+initSmallEM <- function(k, X, Y, fast)
 {
-	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))
-	Gam = matrix(0, n, k)
-	piInit1 = matrix(0,20,k)
-	gamInit1 = array(0, dim=c(n,k,20))
-	LLFinit1 = list()
+  n <- nrow(Y)
+  m <- ncol(Y)
+  p <- ncol(X)
+  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, 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)
-	{
-	  distance_clus = dist(X)
-	  tree_hier = hclust(distance_clus)
-	  Zinit1[,repet] = cutree(tree_hier, k)
+  # 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) {
-			  betaInit1[,,r,repet] = ginv(crossprod(t(X[Z_indice,]))) %*%
-			    crossprod(t(X[Z_indice,]), Y[Z_indice,])
-			} else {
-			betaInit1[,,r,repet] = ginv(crossprod(X[Z_indice,])) %*%
-				crossprod(X[Z_indice,], Y[Z_indice,])
-			}
-			sigmaInit1[,,r,repet] = diag(m)
-			phiInit1[,,r,repet] = betaInit1[,,r,repet] #/ sigmaInit1[,,r,repet]
-			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)
-			}
-			sumGamI = sum(Gam[i,])
-			gamInit1[i,,repet]= Gam[i,] / sumGamI
-		}
-		
-		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
-		LLFinit1[repet] = LLFEessai[length(LLFEessai)]
-	}
+    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) {
+        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, ])) %*% 
+          crossprod(X[Z_indice, ], Y[Z_indice, ])
+      }
+      sigmaInit1[, , r, repet] <- diag(m)
+      phiInit1[, , r, repet] <- betaInit1[, , r, repet]  #/ sigmaInit1[,,r,repet]
+      rhoInit1[, , r, repet] <- solve(sigmaInit1[, , r, repet])
+      piInit1[repet, r] <- mean(Z == r)
+    }
 
-	b = which.max(LLFinit1)
-	phiInit = phiInit1[,,,b]
-	rhoInit = rhoInit1[,,,b]
-	piInit = piInit1[b,]
-	gamInit = gamInit1[,,b]
+    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)
+      }
+      sumGamI <- sum(Gam[i, ])
+      gamInit1[i, , repet] <- Gam[i, ]/sumGamI
+    }
 
-	return (list(phiInit=phiInit, rhoInit=rhoInit, piInit=piInit, gamInit=gamInit))
+    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)
+    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))
 }