1 EMGLLF = function(phiInit,rhoInit,piInit,gamInit,mini,maxi,gamma,lambda,X,Y,tau){
13 S = array(0, dim=c(p,m,k))
17 Gram2 = array(0, dim=c(p,p,k))
18 ps2 = array(0, dim=c(p,m,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))
29 ps1 = array(0, dim=c(n,m,k))
30 nY21 = array(0, dim=c(n,m,k))
34 while(ite <= mini || (ite<= maxi && (dist>= tau || dist2 >= sqrt(tau)))){
38 #calcul associé à Y et X
41 Y2[,mm,r] = sqrt(gam[,r]) .* Y[,mm]
44 X2[i,,r] = X[i,] .* sqrt(gam[i,r])
47 ps2[,mm,r] = crossprod(X2[,,r],Y2[,mm,r])
51 Gram2[j,s,r] = tcrossprod(X2[,j,r], X2[,s,r])
62 b[r] = sum(sum(abs(phi[,,r])))
65 a = sum(gam*t(log(Pi)))
67 #tant que les props sont negatives
69 pi2AllPositive = FALSE
70 while(pi2AllPositive == FALSE){
71 Pi2 = Pi + 0.1^kk * ((1/n)*gam2 - Pi)
75 pi2AllPositive = false;
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)
89 Pi = (Pi+t*(Pi2-Pi)) / sum(Pi+t*(Pi2-Pi))
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
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]))
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)){
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]
113 phi[j,mm,r] = -(n*lambda*(Pi[r]^gamma)+S[j,mm,r])/Gram2[j,j,r]
125 #precompute dot products to numerically adjust their values
126 dotProducts = rep(0,k)
128 dotProducts[r] = tcrossprod(Y[i,]%*%rho[,,r]-X[i,]%*%phi[,,r])
130 shift = 0.5*min(dotProducts)
132 #compute Gam(:,:) using shift determined above
135 Gam[i,r] = Pi[r]*det(rho[,,r])*exp(-0.5*dotProducts[r] + shift)
136 sumLLF1 = sumLLF1 + Gam[i,r]/(2*pi)^(m/2)
138 sumLogLLF2 = sumLogLLF2 + log(sumLLF1)
139 sumGamI = sum(Gam[i,])
141 gam[i,] = Gam[i,] / sumGamI
149 sumPen = sumPen + Pi[r].^gamma^b[r]
151 LLF[ite] = -(1/n)*sumLogLLF2 + lambda*sumPen
156 dist = (LLF[ite]-LLF[ite-1])/(1+abs(LLF[ite]))
158 Dist1=max(max(max((abs(phi-Phi))./(1+abs(phi)))))
159 Dist2=max(max(max((abs(rho-Rho))./(1+abs(rho)))))
160 Dist3=max(max((abs(Pi-PI))./(1+abs(PI))))
161 dist2=max([Dist1,Dist2,Dist3])
167 return(list(phi=phi, rho=rho, Pi=Pi, LLF=LLF, S=S))