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(n,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] ##bon calcul ? idem pour X2 ??...
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] = crossprod(X2[,j,r], X2[,s,r])
62 b[r] = sum(sum(abs(phi[,,r])))
65 a = sum(gam%*%(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)%*%t(b)))<(-1/n*gam2%*%t(log(Pi2))+lambda*(Pi2^gamma)%*%t(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] * (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]))
111 S[j,mm,r] = -rho[mm,mm,r]*ps2[j,mm,r] + phi[v1,mm,r]%*%(Gram2[j,v1,r]) + phi[v2,mm,r]%*%(Gram2[j,v2,r]) #erreur indice
112 if(abs(S[j,mm,r]) <= n*lambda*(Pi[r]^gamma)){
115 if(S[j,mm,r]> n*lambda*(Pi[r]^gamma)){
116 phi[j,mm,r] = (n*lambda*(Pi[r]^gamma)-S[j,mm,r])/Gram2[j,j,r]
118 phi[j,mm,r] = -(n*lambda*(Pi[r]^gamma)+S[j,mm,r])/Gram2[j,j,r]
130 #precompute dot products to numerically adjust their values
131 dotProducts = rep(0,k)
133 dotProducts[r] = tcrossprod(Y[i,]%*%rho[,,r]-X[i,]%*%phi[,,r])
135 shift = 0.5*min(dotProducts)
137 #compute Gam(:,:) using shift determined above
140 Gam[i,r] = Pi[r]*det(rho[,,r])*exp(-0.5*dotProducts[r] + shift)
141 sumLLF1 = sumLLF1 + Gam[i,r]/(2*pi)^(m/2)
143 sumLogLLF2 = sumLogLLF2 + log(sumLLF1)
144 sumGamI = sum(Gam[i,])
146 gam[i,] = Gam[i,] / sumGamI
154 sumPen = sumPen + Pi[r]^gamma^b[r]
156 LLF[ite] = -(1/n)*sumLogLLF2 + lambda*sumPen
161 dist = (LLF[ite]-LLF[ite-1])/(1+abs(LLF[ite]))
163 Dist1=max(max(max((abs(phi-Phi))/(1+abs(phi)))))
164 Dist2=max(max(max((abs(rho-Rho))/(1+abs(rho)))))
165 Dist3=max(max((abs(Pi-PI))/(1+abs(PI))))
166 dist2=max(c(Dist1,Dist2,Dist3))
172 return(list("phi"=phi, "rho"=rho, "pi"=Pi, "LLF"=LLF, "S"=S))