test/generate_test_data/EMGLLF.R

   1 EMGLLF_R = function(phiInit,rhoInit,piInit,gamInit,mini,maxi,gamma,lambda,X,Y,tau)
   2 {
   3     # Matrix dimensions
   4     n = dim(X)[1]
   5     p = dim(phiInit)[1]
   6     m = dim(phiInit)[2]
   7     k = dim(phiInit)[3]
   8
   9     # Outputs
  10     phi = phiInit
  11     rho = rhoInit
  12     pi = piInit
  13     llh = -Inf
  14     S = array(0, dim=c(p,m,k))
  15
  16     # Algorithm variables
  17     gam = gamInit
  18     Gram2 = array(0, dim=c(p,p,k))
  19     ps2 = array(0, dim=c(p,m,k))
  20     X2 = array(0, dim=c(n,p,k))
  21     Y2 = array(0, dim=c(n,m,k))
  22     EPS = 1e-15
  23
  24     for (ite in 1:maxi)
  25     {
  26         # Remember last pi,rho,phi values for exit condition in the end of loop
  27         Phi = phi
  28         Rho = rho
  29         Pi = pi
  30
  31         # Calcul associé à Y et X
  32         for (r in 1:k)
  33         {
  34             for (mm in 1:m)
  35                 Y2[,mm,r] = sqrt(gam[,r]) * Y[,mm]
  36             for (i in 1:n)
  37                 X2[i,,r] = sqrt(gam[i,r]) * X[i,]
  38             for (mm in 1:m)
  39                 ps2[,mm,r] = crossprod(X2[,,r],Y2[,mm,r])
  40             for (j in 1:p)
  41             {
  42                 for (s in 1:p)
  43                     Gram2[j,s,r] = crossprod(X2[,j,r], X2[,s,r])
  44             }
  45         }
  46
  47         ##########
  48         #Etape M #
  49         ##########
  50
  51         # Pour pi
  52         b = sapply( 1:k, function(r) sum(abs(phi[,,r])) )
  53         gam2 = colSums(gam)
  54         a = sum(gam %*% log(pi))
  55
  56         # Tant que les props sont negatives
  57         kk = 0
  58         pi2AllPositive = FALSE
  59         while (!pi2AllPositive)
  60         {
  61             pi2 = pi + 0.1^kk * ((1/n)*gam2 - pi)
  62             pi2AllPositive = all(pi2 >= 0)
  63             kk = kk+1
  64         }
  65
  66         # t(m) la plus grande valeur dans la grille O.1^k tel que ce soit décroissante ou constante
  67         while( kk < 1000 && -a/n + lambda * sum(pi^gamma * b) <
  68             -sum(gam2 * log(pi2))/n + lambda * sum(pi2^gamma * b) )
  69         {
  70             pi2 = pi + 0.1^kk * (1/n*gam2 - pi)
  71             kk = kk + 1
  72         }
  73         t = 0.1^kk
  74         pi = (pi + t*(pi2-pi)) / sum(pi + t*(pi2-pi))
  75
  76         #Pour phi et rho
  77         for (r in 1:k)
  78         {
  79             for (mm in 1:m)
  80             {
  81                 ps = 0
  82                 for (i in 1:n)
  83                     ps = ps + Y2[i,mm,r] * sum(X2[i,,r] * phi[,mm,r])
  84                 nY2 = sum(Y2[,mm,r]^2)
  85                 rho[mm,mm,r] = (ps+sqrt(ps^2+4*nY2*gam2[r])) / (2*nY2)
  86             }
  87         }
  88
  89         for (r in 1:k)
  90         {
  91             for (j in 1:p)
  92             {
  93                 for (mm in 1:m)
  94                 {
  95                     S[j,mm,r] = -rho[mm,mm,r]*ps2[j,mm,r] + sum(phi[-j,mm,r] * Gram2[j,-j,r])
  96                     if (abs(S[j,mm,r]) <= n*lambda*(pi[r]^gamma))
  97                         phi[j,mm,r]=0
  98                     else if(S[j,mm,r] > n*lambda*(pi[r]^gamma))
  99                         phi[j,mm,r] = (n*lambda*(pi[r]^gamma)-S[j,mm,r]) / Gram2[j,j,r]
 100                     else
 101                         phi[j,mm,r] = -(n*lambda*(pi[r]^gamma)+S[j,mm,r]) / Gram2[j,j,r]
 102                 }
 103             }
 104         }
 105
 106         ##########
 107         #Etape E #
 108         ##########
 109
 110         # Precompute det(rho[,,r]) for r in 1...k
 111         detRho = sapply(1:k, function(r) det(rho[,,r]))
 112
 113         sumLogLLH = 0
 114         for (i in 1:n)
 115         {
 116             # Update gam[,]
 117             sumGamI = 0
 118             for (r in 1:k)
 119             {
 120                 gam[i,r] = pi[r]*exp(-0.5*sum((Y[i,]%*%rho[,,r]-X[i,]%*%phi[,,r])^2))*detRho[r]
 121                 sumGamI = sumGamI + gam[i,r]
 122             }
 123             sumLogLLH = sumLogLLH + log(sumGamI) - log((2*base::pi)^(m/2))
 124             if (sumGamI > EPS) #else: gam[i,] is already ~=0
 125                 gam[i,] = gam[i,] / sumGamI
 126         }
 127
 128         sumPen = sum(pi^gamma * b)
 129         last_llh = llh
 130         llh = -sumLogLLH/n + lambda*sumPen
 131         dist = ifelse( ite == 1, llh, (llh-last_llh) / (1+abs(llh)) )
 132         Dist1 = max( (abs(phi-Phi)) / (1+abs(phi)) )
 133         Dist2 = max( (abs(rho-Rho)) / (1+abs(rho)) )
 134         Dist3 = max( (abs(pi-Pi)) / (1+abs(Pi)) )
 135         dist2 = max(Dist1,Dist2,Dist3)
 136
 137         if (ite >= mini && (dist >= tau || dist2 >= sqrt(tau)))
 138             break
 139     }
 140
 141     affec = apply(gam, 1, which.max)
 142     list( "phi"=phi, "rho"=rho, "pi"=pi, "llh"=llh, "S"=S, "affec"=affec )
 143 }