src/test/generate_test_data/helpers/EMGLLF.R

   1 EMGLLF = function(phiInit,rhoInit,piInit,gamInit,mini,maxi,gamma,lambda,X,Y,tau){
   2   #matrix dimensions
   3   n = dim(X)[1]
   4   p = dim[phiInit][1]
   5   m = dim[phiInit][2]
   6   k = dim[phiInit][3]
   7
   8   #init outputs
   9   phi = phiInit
  10   rho = rhoInit
  11   Pi = piInit
  12   LLF = rep(0, maxi)
  13   S = array(0, dim=c(p,m,k))
  14
  15
  16   gam = gamInit
  17   Gram2 = array(0, dim=c(p,p,k))
  18   ps2 = array(0, dim=c(p,m,k))
  19   b = rep(0, k)
  20   pen = matrix(0, maxi, k)
  21   X2 = array(0, dim=c(n,p,k))
  22   Y2 = array(0, dim=c(p,m,k))
  23   dist = 0
  24   dist2 = 0
  25   ite = 1
  26   Pi2 = rep(0, k)
  27   ps = matrix(0, m,k)
  28   nY2 = matrix(0, m,k)
  29   ps1 = array(0, dim=c(n,m,k))
  30   nY21 = array(0, dim=c(n,m,k))
  31   Gam = matrix(0, n,k)
  32   EPS = 1E-15
  33
  34   while(ite <= mini || (ite<= maxi && (dist>= tau || dist2 >= sqrt(tau)))){
  35     Phi = phi
  36     Rho = rho
  37     PI = Pi
  38     #calcul associé à Y et X
  39     for(r in 1:k){
  40       for(mm in 1:m){
  41         Y2[,mm,r] = sqrt(gam[,r]) .* Y[,mm]
  42       }
  43       for(i in 1:n){
  44         X2[i,,r] = X[i,] .* sqrt(gam[i,r])
  45       }
  46       for(mm in 1:m){
  47         ps2[,mm,r] = crossprod(X2[,,r],Y2[,mm,r])
  48       }
  49       for(j in 1:p){
  50         for(s in 1:p){
  51           Gram2[j,s,r] = tcrossprod(X2[,j,r], X2[,s,r])
  52         }
  53       }
  54     }
  55
  56     ##########
  57     #Etape M #
  58     ##########
  59
  60     #pour pi
  61     for(r in 1:k){
  62       b[r] = sum(sum(abs(phi[,,r])))
  63     }
  64     gam2 = sum(gam[1,])  #BIG DOUTE
  65     a = sum(gam*t(log(Pi)))
  66
  67     #tant que les props sont negatives
  68     kk = 0
  69     pi2AllPositive = FALSE
  70     while(pi2AllPositive == FALSE){
  71       pi2 = pi + 0.1^kk * ((1/n)*gam2 - pi)
  72       pi2AllPositive = TRUE
  73       for(r in 1:k){
  74         if(pi2[r] < 0){
  75           pi2AllPositive = false;
  76           break
  77         }
  78       }
  79       kk = kk+1
  80     }
  81
  82     #t[m]la plus grande valeur dans la grille O.1^k tel que ce soit
  83     #décroissante ou constante
  84     while((-1/n*a+lambda*((pi.^gamma)*b))<(-1/n*gam2*t(log(pi2))+lambda.*(pi2.^gamma)*b) && kk<1000){
  85       pi2 = pi+0.1^kk*(1/n*gam2-pi)
  86       kk = kk+1
  87     }
  88     t = 0.1^(kk)
  89     pi = (pi+t*(pi2-pi)) / sum(pi+t*(pi2-pi))
  90
  91     #Pour phi et rho
  92     for(r in 1:k){
  93       for(mm in 1:m){
  94         for(i in 1:n){
  95           ps1[i,mm,r] = Y2[i,mm,r] * dot(X2(i,:,r), phi(:,mm,r))
  96           nY21[i,mm,r] = (Y2[i,mm,r])^2
  97         }
  98         ps[mm,r] = sum(ps1(:,mm,r));
  99         nY2[mm,r] = sum(nY21(:,mm,r));
 100         rho[mm,mm,r] = ((ps[mm,r]+sqrt(ps[mm,r]^2+4*nY2[mm,r]*(gam2[r])))/(2*nY2[mm,r]))
 101       }
 102     }
 103     for(r in 1:k){
 104       for(j in 1:p){
 105         for(mm in 1:m){
 106           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])
 107           if(abs(S(j,mm,r)) <= n*lambda*(pi(r)^gamma))
 108             phi[j,mm,r]=0
 109           else{
 110             if(S[j,mm,r]> n*lambda*(Pi[r]^gamma))
 111               phi[j,mm,r] = (n*lambda*(Pi[r]^gamma)-S[j,mm,r])/Gram2[j,j,r]
 112           else
 113             phi[j,mm,r] = -(n*lambda*(Pi[r]^gamma)+S[j,mm,r])/Gram2[j,j,r]
 114           }
 115         }
 116       }
 117     }
 118
 119     ##########
 120     #Etape E #
 121     ##########
 122     sumLogLLF2 = 0
 123     for(i in 1:n){
 124       #precompute dot products to numerically adjust their values
 125       dotProducts = rep(0,k)
 126       for(r in 1:k){
 127         dotProducts[r] = tcrossprod(Y[i,]%*%rho[,,r]-X[i,]%*%phi[,,r])
 128       }
 129       shift = 0.5*min(dotProducts)
 130
 131       #compute Gam(:,:) using shift determined above
 132       sumLLF1 = 0.0;
 133       for(r in 1:k){
 134         Gam[i,r] = Pi[r]*det(rho[,,r])*exp(-0.5*dotProducts[r] + shift)
 135         sumLLF1 = sumLLF1 + Gam[i,r]/(2*pi)^(m/2)
 136       }
 137       sumLogLLF2 = sumLogLLF2 + log(sumLLF1)
 138       sumGamI = sum(Gam[i,])
 139       if(sumGamI > EPS)
 140         gam[i,] = Gam[i,] / sumGamI
 141       else
 142         gam[i,] = rep(0,k)
 143     }
 144
 145
 146     sumPen = 0
 147     for(r in 1:k){
 148       sumPen = sumPen + Pi[r].^gamma^b[r]
 149     }
 150     LLF[ite] = -(1/n)*sumLogLLF2 + lambda*sumPen
 151
 152     if(ite == 1)
 153       dist = LLF[ite]
 154     else
 155       dist = (LLF[ite]-LLF[ite-1])/(1+abs(LLF[ite]))
 156
 157     Dist1=max(max(max((abs(phi-Phi))./(1+abs(phi)))))
 158     Dist2=max(max(max((abs(rho-Rho))./(1+abs(rho)))))
 159     Dist3=max(max((abs(Pi-PI))./(1+abs(PI))))
 160     dist2=max([Dist1,Dist2,Dist3])
 161
 162     ite=ite+1
 163   }
 164
 165   Pi = transpose(Pi)
 166   return(list(phi=phi, rho=rho, Pi=Pi, LLF=LLF, S=S))
 167 }