EMGrank + typo

diff --git a/src/test/generate_test_data/helpers/EMGLLF.R b/src/test/generate_test_data/helpers/EMGLLF.R
index f108a38..94a917d 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]
 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
   
   #init outputs
   phi = phiInit


@@ -68,10 +68,10 @@ EMGLLF = function(phiInit,rhoInit,piInit,gamInit,mini,maxi,gamma,lambda,X,Y,tau)
     kk = 0
     pi2AllPositive = FALSE
     while(pi2AllPositive == FALSE){
     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){
       pi2AllPositive = TRUE
       for(r in 1:k){

-        if(pi2[r] < 0){
+        if(Pi2[r] < 0){

           pi2AllPositive = false;
           break
         }
           pi2AllPositive = false;
           break
         }


@@ -81,12 +81,12 @@ 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
     
     #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)*b))<(-1/n*gam2*t(log(Pi2))+lambda.*(Pi2.^gamma)*b) && kk<1000){
+      Pi2 = Pi+0.1^kk*(1/n*gam2-Pi)

       kk = kk+1
     }
     t = 0.1^(kk)
       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){
     
     #Pour phi et rho
     for(r in 1:k){


@@ -104,13 +104,14 @@ EMGLLF = function(phiInit,rhoInit,piInit,gamInit,mini,maxi,gamma,lambda,X,Y,tau)
       for(j in 1:p){
         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])
       for(j in 1:p){
         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))
+          if(abs(S(j,mm,r)) <= n*lambda*(Pi[r]^gamma)){

             phi[j,mm,r]=0
             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]
               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]
+            }

           }
         }
       }
           }
         }
       }

diff --git a/src/test/generate_test_data/helpers/EMGrank.R b/src/test/generate_test_data/helpers/EMGrank.R
new file mode 100644 (file)
index 0000000..0ee2d1f
--- /dev/null
+++ b/src/test/generate_test_data/helpers/EMGrank.R
@@ -0,0 +1,72 @@
+EMGLLF = function(Pi, Rho, mini, maxi, X, Y, tau, rank){
+  #matrix dimensions
+  n = dim(X)[1]
+  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)
+  Pi = piInit
+  LLF = 0
+  
+  #local variables
+  Phi = array(0, dim=c(p,m,k))
+  deltaPhi = c(0)
+  sumDeltaPhi = 0
+  deltaPhiBufferSize = 20
+  
+  #main loop
+  ite = 1
+  while(ite<=mini || (ite<=maxi && sumDeltaPhi>tau)){
+    #M step: Mise à jour de Beta (et donc phi)
+    for(r in 1:k){
+      Z_bin = vec_bin(Z,r)
+      Z_vec = Z_bin$vec #vecteur 0 et 1 aux endroits o? Z==r
+      Z_indice = Z_bin$indice 
+      if(sum(Z_indice) == 0){
+        next
+      }
+      #U,S,V = SVD of (t(Xr)Xr)^{-1} * t(Xr) * Yr
+      [U,S,V] = svd(ginv(crossprod(X[Z_indice,]))%*% (X[Z_indice,])%*%Y[Z_indice,]      )
+      #Set m-rank(r) singular values to zero, and recompose
+      #best rank(r) approximation of the initial product
+      S[rank(r)+1:end,] = 0
+      phi[,,r] = U %*%S%*%t(V)%*%Rho[,,r]
+    }
+  
+  #Etape E et calcul de LLF
+  sumLogLLF2 = 0
+  for(i in 1:n){
+    sumLLF1 = 0
+    maxLogGamIR = -Inf
+    for(r in 1:k){
+      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){
+        Z[i] = r
+        maxLogGamIR = logGamIR
+      }
+    sumLLF1 = sumLLF1 + exp(logGamIR) / (2*pi)^(m/2)
+    }
+    sumLogLLF2 = sumLogLLF2 + log(sumLLF1)
+  }
+  
+  LLF = -1/n * sumLogLLF2
+  
+  #update distance parameter to check algorithm convergence (delta(phi, Phi))
+  deltaPhi = c(deltaPhi, max(max(max((abs(phi-Phi))/(1+abs(phi))))) )
+  if(length(deltaPhi) > deltaPhiBufferSize){
+    deltaPhi = deltaPhi[2:length(deltaPhi)]
+  }
+  sumDeltaPhi = sum(abs(deltaPhi))
+  
+  #update other local variables
+  Phi = phi
+  ite = ite+1
+  
+  }
+  return(list(phi=phi, LLF=LLF))
+}
\ No newline at end of file

diff --git a/src/test/generate_test_data/helpers/checkOutput.R b/src/test/generate_test_data/helpers/checkOutput.R
index 187536a..ae2fbdf 100644 (file)
--- a/src/test/generate_test_data/helpers/checkOutput.R
+++ b/src/test/generate_test_data/helpers/checkOutput.R
@@ -2,11 +2,9 @@ checkOutput = function(varName, array, refArray, tol)
 {
        print(paste("Checking ",varName,sep=""))
        maxError = max(abs(array - refArray))
 {
        print(paste("Checking ",varName,sep=""))
        maxError = max(abs(array - refArray))

-       if(maxError >= tol)
-       {
+       if(maxError >= tol){

                print(paste("Inaccuracy: max(abs(error)) = ",maxError," >= ",tol,sep=""))
                print(paste("Inaccuracy: max(abs(error)) = ",maxError," >= ",tol,sep=""))

-       } else
-       {
+       } else{

                print("OK")
        }
 }
                print("OK")
        }
 }

diff --git a/src/test/generate_test_data/helpers/covariance.R b/src/test/generate_test_data/helpers/covariance.R
index 15cd693..09a9ec5 100644 (file)
--- a/src/test/generate_test_data/helpers/covariance.R
+++ b/src/test/generate_test_data/helpers/covariance.R
@@ -1,8 +1,8 @@
 covariance = function(p,a)
 {
        A = matrix(a, p,p)
 covariance = function(p,a)
 {
        A = matrix(a, p,p)

-       for (i in 1:p)
+       for (i in 1:p){

                A[i,] = A[i,]^abs(i-(1:p))
                A[i,] = A[i,]^abs(i-(1:p))

-
+       }

        return (A)
 }
        return (A)
 }