+++ /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