--- /dev/null
+function[phi,rho,pi,LLF,S] = EMGLLF(phiInit,rhoInit,piInit,gamInit,mini,maxi,gamma,lambda,X,Y,tau)\r
+\r
+ %Get matrices dimensions\r
+ PI = 4.0 * atan(1.0);\r
+ n = size(X, 1);\r
+ [p,m,k] = size(phiInit);\r
+\r
+ %Initialize outputs\r
+ phi = phiInit;\r
+ rho = rhoInit;\r
+ pi = piInit;\r
+ LLF = zeros(maxi,1);\r
+ S = zeros(p,m,k);\r
+\r
+ %Other local variables\r
+ %NOTE: variables order is always n,p,m,k\r
+ gam = gamInit;\r
+ Gram2 = zeros(p,p,k);\r
+ ps2 = zeros(p,m,k);\r
+ b = zeros(k,1);\r
+ pen = zeros(maxi,k);\r
+ X2 = zeros(n,p,k);\r
+ Y2 = zeros(n,m,k);\r
+ dist = 0;\r
+ dist2 = 0;\r
+ ite = 1;\r
+ pi2 = zeros(k,1);\r
+ ps = zeros(m,k);\r
+ nY2 = zeros(m,k);\r
+ ps1 = zeros(n,m,k);\r
+ nY21 = zeros(n,m,k);\r
+ Gam = zeros(n,k);\r
+ EPS = 1e-15;\r
+\r
+ while ite<=mini || (ite<=maxi && (dist>=tau || dist2>=sqrt(tau)))\r
+\r
+ Phi = phi;\r
+ Rho = rho;\r
+ Pi = pi;\r
+\r
+ %Calculs associés à Y et X\r
+ for r=1:k\r
+ for mm=1:m\r
+ Y2(:,mm,r) = sqrt(gam(:,r)) .* Y(:,mm);\r
+ end\r
+ for i=1:n\r
+ X2(i,:,r) = X(i,:) .* sqrt(gam(i,r));\r
+ end\r
+ for mm=1:m\r
+ ps2(:,mm,r) = transpose(X2(:,:,r)) * Y2(:,mm,r);\r
+ end\r
+ for j=1:p\r
+ for s=1:p\r
+ Gram2(j,s,r) = dot(X2(:,j,r), X2(:,s,r));\r
+ end\r
+ end\r
+ end\r
+\r
+ %%%%%%%%%%\r
+ %Etape M %\r
+ %%%%%%%%%%\r
+\r
+ %Pour pi\r
+ for r=1:k\r
+ b(r) = sum(sum(abs(phi(:,:,r))));\r
+ end\r
+ gam2 = sum(gam,1);\r
+ a = sum(gam*transpose(log(pi)));\r
+\r
+ %tant que les proportions sont negatives\r
+ kk = 0;\r
+ pi2AllPositive = false;\r
+ while ~pi2AllPositive\r
+ pi2 = pi + 0.1^kk * ((1/n)*gam2 - pi);\r
+ pi2AllPositive = true;\r
+ for r=1:k\r
+ if pi2(r) < 0\r
+ pi2AllPositive = false;\r
+ break;\r
+ end\r
+ end\r
+ kk = kk+1;\r
+ end\r
+\r
+ %t(m) la plus grande valeur dans la grille O.1^k tel que ce soit\r
+ %décroissante ou constante\r
+ while (-1/n*a+lambda*((pi.^gamma)*b))<(-1/n*gam2*transpose(log(pi2))+lambda.*(pi2.^gamma)*b) && kk<1000\r
+ pi2 = pi+0.1^kk*(1/n*gam2-pi);\r
+ kk = kk+1;\r
+ end\r
+ t = 0.1^(kk);\r
+ pi = (pi+t*(pi2-pi)) / sum(pi+t*(pi2-pi));\r
+\r
+ %Pour phi et rho\r
+ for r=1:k\r
+ for mm=1:m\r
+ for i=1:n\r
+ ps1(i,mm,r) = Y2(i,mm,r) * dot(X2(i,:,r), phi(:,mm,r));\r
+ nY21(i,mm,r) = (Y2(i,mm,r))^2;\r
+ end\r
+ ps(mm,r) = sum(ps1(:,mm,r));\r
+ nY2(mm,r) = sum(nY21(:,mm,r));\r
+ rho(mm,mm,r) = ((ps(mm,r)+sqrt(ps(mm,r)^2+4*nY2(mm,r)*(gam2(r))))/(2*nY2(mm,r)));\r
+ end\r
+ end\r
+ for r=1:k\r
+ for j=1:p\r
+ for mm=1:m\r
+ 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)')...\r
+ + dot(phi(j+1:p,mm,r),Gram2(j,j+1:p,r)');\r
+ if abs(S(j,mm,r)) <= n*lambda*(pi(r)^gamma)\r
+ phi(j,mm,r)=0;\r
+ else\r
+ if S(j,mm,r)> n*lambda*(pi(r)^gamma)\r
+ phi(j,mm,r)=(n*lambda*(pi(r)^gamma)-S(j,mm,r))/Gram2(j,j,r);\r
+ else\r
+ phi(j,mm,r)=-(n*lambda*(pi(r)^gamma)+S(j,mm,r))/Gram2(j,j,r);\r
+ end\r
+ end\r
+ end\r
+ end\r
+ end\r
+\r
+ %%%%%%%%%%\r
+ %Etape E %\r
+ %%%%%%%%%%\r
+\r
+ sumLogLLF2 = 0.0;\r
+ for i=1:n\r
+ %precompute dot products to numerically adjust their values\r
+ dotProducts = zeros(k,1);\r
+ for r=1:k\r
+ dotProducts(r)= (Y(i,:)*rho(:,:,r)-X(i,:)*phi(:,:,r)) * transpose(Y(i,:)*rho(:,:,r)-X(i,:)*phi(:,:,r));\r
+ end\r
+ shift = 0.5*min(dotProducts);\r
+\r
+ %compute Gam(:,:) using shift determined above\r
+ sumLLF1 = 0.0;\r
+ for r=1:k\r
+ Gam(i,r) = pi(r)*det(rho(:,:,r))*exp(-0.5*dotProducts(r) + shift);\r
+ sumLLF1 = sumLLF1 + Gam(i,r)/(2*PI)^(m/2);\r
+ end\r
+ sumLogLLF2 = sumLogLLF2 + log(sumLLF1);\r
+ sumGamI = sum(Gam(i,:));\r
+ if sumGamI > EPS\r
+ gam(i,:) = Gam(i,:) / sumGamI;\r
+ else\r
+ gam(i,:) = zeros(k,1);\r
+ end\r
+ end\r
+\r
+ sumPen = 0.0;\r
+ for r=1:k\r
+ sumPen = sumPen + pi(r).^gamma .* b(r);\r
+ end\r
+ LLF(ite) = -(1/n)*sumLogLLF2 + lambda*sumPen;\r
+\r
+ if ite == 1\r
+ dist = LLF(ite);\r
+ else\r
+ dist = (LLF(ite)-LLF(ite-1))/(1+abs(LLF(ite)));\r
+ end\r
+\r
+ Dist1=max(max(max((abs(phi-Phi))./(1+abs(phi)))));\r
+ Dist2=max(max(max((abs(rho-Rho))./(1+abs(rho)))));\r
+ Dist3=max(max((abs(pi-Pi))./(1+abs(Pi))));\r
+ dist2=max([Dist1,Dist2,Dist3]);\r
+\r
+ ite=ite+1;\r
+ end\r
+\r
+ pi = transpose(pi);\r
+\r
+end\r
projects / valse.git / blobdiff