| 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 | } |