correction emgllf.R

[valse.git] / src / test / generate_test_data / helpers / EMGLLF.R

diff --git a/src/test/generate_test_data/helpers/EMGLLF.R b/src/test/generate_test_data/helpers/EMGLLF.R
index f108a38..b0d95a3 100644 (file)
--- a/src/test/generate_test_data/helpers/EMGLLF.R
+++ b/src/test/generate_test_data/helpers/EMGLLF.R
@@ -1,9 +1,9 @@
 EMGLLF = function(phiInit,rhoInit,piInit,gamInit,mini,maxi,gamma,lambda,X,Y,tau){
   #matrix dimensions
   n = dim(X)[1]
-  p = dim[phiInit][1]
-  m = dim[phiInit][2]
-  k = dim[phiInit][3]
+  p = dim(phiInit)[1]
+  m = dim(phiInit)[2]
+  k = dim(phiInit)[3]
   
   #init outputs
   phi = phiInit
@@ -19,7 +19,7 @@ EMGLLF = function(phiInit,rhoInit,piInit,gamInit,mini,maxi,gamma,lambda,X,Y,tau)
   b = rep(0, k)
   pen = matrix(0, maxi, k)
   X2 = array(0, dim=c(n,p,k))
-  Y2 = array(0, dim=c(p,m,k))
+  Y2 = array(0, dim=c(n,m,k))
   dist = 0
   dist2 = 0
   ite = 1
@@ -38,17 +38,17 @@ EMGLLF = function(phiInit,rhoInit,piInit,gamInit,mini,maxi,gamma,lambda,X,Y,tau)
     #calcul associé à Y et X
     for(r in 1:k){
       for(mm in 1:m){
-        Y2[,mm,r] = sqrt(gam[,r]) .* Y[,mm]
+        Y2[,mm,r] = sqrt(gam[,r]) * Y[,mm]  ##bon calcul ? idem pour X2 ??...
       }
       for(i in 1:n){
-        X2[i,,r] = X[i,] .* sqrt(gam[i,r])
+        X2[i,,r] = X[i,] *sqrt(gam[i,r])
       }
       for(mm in 1:m){
         ps2[,mm,r] = crossprod(X2[,,r],Y2[,mm,r])
       }
       for(j in 1:p){
         for(s in 1:p){
-          Gram2[j,s,r] = tcrossprod(X2[,j,r], X2[,s,r])
+          Gram2[j,s,r] = crossprod(X2[,j,r], X2[,s,r])
         }
       }
     }
@@ -61,17 +61,17 @@ EMGLLF = function(phiInit,rhoInit,piInit,gamInit,mini,maxi,gamma,lambda,X,Y,tau)
     for(r in 1:k){
       b[r] = sum(sum(abs(phi[,,r])))
     }
-    gam2 = sum(gam[1,])  #BIG DOUTE
-    a = sum(gam*t(log(Pi)))
+    gam2 = colSums(gam)
+    a = sum(gam%*%(log(Pi)))
     
     #tant que les props sont negatives
     kk = 0
     pi2AllPositive = FALSE
     while(pi2AllPositive == FALSE){
-      pi2 = pi + 0.1^kk * ((1/n)*gam2 - pi)
+      Pi2 = Pi + 0.1^kk * ((1/n)*gam2 - Pi)
       pi2AllPositive = TRUE
       for(r in 1:k){
-        if(pi2[r] < 0){
+        if(Pi2[r] < 0){
           pi2AllPositive = false;
           break
         }
@@ -81,36 +81,42 @@ EMGLLF = function(phiInit,rhoInit,piInit,gamInit,mini,maxi,gamma,lambda,X,Y,tau)
     
     #t[m]la plus grande valeur dans la grille O.1^k tel que ce soit
     #décroissante ou constante
-    while((-1/n*a+lambda*((pi.^gamma)*b))<(-1/n*gam2*t(log(pi2))+lambda.*(pi2.^gamma)*b) && kk<1000){
-      pi2 = pi+0.1^kk*(1/n*gam2-pi)
+    while((-1/n*a+lambda*((Pi^gamma)%*%t(b)))<(-1/n*gam2%*%t(log(Pi2))+lambda*(Pi2^gamma)%*%t(b)) && kk<1000){
+      Pi2 = Pi+0.1^kk*(1/n*gam2-Pi)
       kk = kk+1
     }
     t = 0.1^(kk)
-    pi = (pi+t*(pi2-pi)) / sum(pi+t*(pi2-pi))
+    Pi = (Pi+t*(Pi2-Pi)) / sum(Pi+t*(Pi2-Pi))
     
     #Pour phi et rho
     for(r in 1:k){
       for(mm in 1:m){
         for(i in 1:n){
-          ps1[i,mm,r] = Y2[i,mm,r] * dot(X2(i,:,r), phi(:,mm,r))
+          ps1[i,mm,r] = Y2[i,mm,r] * (X2[i,,r]%*%(phi[,mm,r]))
           nY21[i,mm,r] = (Y2[i,mm,r])^2
         }
-        ps[mm,r] = sum(ps1(:,mm,r));
-        nY2[mm,r] = sum(nY21(:,mm,r));
+        ps[mm,r] = sum(ps1[,mm,r])
+        nY2[mm,r] = sum(nY21[,mm,r])
         rho[mm,mm,r] = ((ps[mm,r]+sqrt(ps[mm,r]^2+4*nY2[mm,r]*(gam2[r])))/(2*nY2[mm,r]))
       }
     }
     for(r in 1:k){
-      for(j in 1:p){
+      p1 = p-1
+      for(j in 1:p1){
         for(mm in 1:m){
-          S[j,mm,r] = -rho[mm,mm,r]*ps2[j,mm,r] + dot(phi[1:j-1,mm,r],Gram2[j,1:j-1,r])  + dot(phi[j+1:p,mm,r],Gram2[j,j+1:p,r])
-          if(abs(S(j,mm,r)) <= n*lambda*(pi(r)^gamma))
+          j1 = j-1
+          j2 = j+1
+          v1 = c(1:j1)
+          v2 = c(j2:p)
+          S[j,mm,r] = -rho[mm,mm,r]*ps2[j,mm,r] + phi[v1,mm,r]%*%(Gram2[j,v1,r])  + phi[v2,mm,r]%*%(Gram2[j,v2,r])  #erreur indice
+          if(abs(S[j,mm,r]) <= n*lambda*(Pi[r]^gamma)){
             phi[j,mm,r]=0
-          else{
-            if(S[j,mm,r]> n*lambda*(Pi[r]^gamma))
+          }else{
+            if(S[j,mm,r]> n*lambda*(Pi[r]^gamma)){
               phi[j,mm,r] = (n*lambda*(Pi[r]^gamma)-S[j,mm,r])/Gram2[j,j,r]
-          else
-            phi[j,mm,r] = -(n*lambda*(Pi[r]^gamma)+S[j,mm,r])/Gram2[j,j,r]
+            }else{
+              phi[j,mm,r] = -(n*lambda*(Pi[r]^gamma)+S[j,mm,r])/Gram2[j,j,r]
+            }
           }
         }
       }
@@ -145,7 +151,7 @@ EMGLLF = function(phiInit,rhoInit,piInit,gamInit,mini,maxi,gamma,lambda,X,Y,tau)
     
     sumPen = 0
     for(r in 1:k){
-      sumPen = sumPen + Pi[r].^gamma^b[r]
+      sumPen = sumPen + Pi[r]^gamma^b[r]
     }
     LLF[ite] = -(1/n)*sumLogLLF2 + lambda*sumPen
     
@@ -154,14 +160,14 @@ EMGLLF = function(phiInit,rhoInit,piInit,gamInit,mini,maxi,gamma,lambda,X,Y,tau)
     else
       dist = (LLF[ite]-LLF[ite-1])/(1+abs(LLF[ite]))
     
-    Dist1=max(max(max((abs(phi-Phi))./(1+abs(phi)))))
-    Dist2=max(max(max((abs(rho-Rho))./(1+abs(rho)))))
-    Dist3=max(max((abs(Pi-PI))./(1+abs(PI))))
-    dist2=max([Dist1,Dist2,Dist3])
+    Dist1=max(max(max((abs(phi-Phi))/(1+abs(phi)))))
+    Dist2=max(max(max((abs(rho-Rho))/(1+abs(rho)))))
+    Dist3=max(max((abs(Pi-PI))/(1+abs(PI))))
+    dist2=max(c(Dist1,Dist2,Dist3))
     
     ite=ite+1
   }
     
-  Pi = transpose(Pi)
+  Pi = t(Pi)
   return(list(phi=phi, rho=rho, Pi=Pi, LLF=LLF, S=S))
-}
\ No newline at end of file
+}