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