--- /dev/null
+function[phi,LLF] = EMGrank(Pi,Rho,mini,maxi,X,Y,tau,rank)\r
+\r
+ % get matrix sizes\r
+ [~,m,k] = size(Rho);\r
+ [n,p] = size(X);\r
+\r
+ % allocate output matrices\r
+ phi = zeros(p,m,k);\r
+ Z = ones(n,1,'int64');\r
+ LLF = 0.0;\r
+\r
+ % local variables\r
+ Phi = zeros(p,m,k);\r
+ deltaPhi = 0.0;\r
+ deltaPhi = [];\r
+ sumDeltaPhi = 0.0;\r
+ deltaPhiBufferSize = 20;\r
+\r
+ %main loop (at least mini iterations)\r
+ ite = int64(1);\r
+ while ite<=mini || (ite<=maxi && sumDeltaPhi>tau)\r
+\r
+ %M step: Mise à jour de Beta (et donc phi)\r
+ for r=1:k\r
+ if (sum(Z==r) == 0)\r
+ continue;\r
+ end\r
+ %U,S,V = SVD of (t(Xr)Xr)^{-1} * t(Xr) * Yr\r
+ [U,S,V] = svd(pinv(transpose(X(Z==r,:))*X(Z==r,:))*transpose(X(Z==r,:))*Y(Z==r,:));\r
+ %Set m-rank(r) singular values to zero, and recompose\r
+ %best rank(r) approximation of the initial product\r
+ S(rank(r)+1:end,:) = 0;\r
+ phi(:,:,r) = U * S * transpose(V) * Rho(:,:,r);\r
+ end\r
+\r
+ %Etape E et calcul de LLF\r
+ sumLogLLF2 = 0.0;\r
+ for i=1:n\r
+ sumLLF1 = 0.0;\r
+ maxLogGamIR = -Inf;\r
+ for r=1:k\r
+ dotProduct = (Y(i,:)*Rho(:,:,r)-X(i,:)*phi(:,:,r)) * transpose(Y(i,:)*Rho(:,:,r)-X(i,:)*phi(:,:,r));\r
+ logGamIR = log(Pi(r)) + log(det(Rho(:,:,r))) - 0.5*dotProduct;\r
+ %Z(i) = index of max (gam(i,:))\r
+ if logGamIR > maxLogGamIR\r
+ Z(i) = r;\r
+ maxLogGamIR = logGamIR;\r
+ end\r
+ sumLLF1 = sumLLF1 + exp(logGamIR) / (2*pi)^(m/2);\r
+ end\r
+ sumLogLLF2 = sumLogLLF2 + log(sumLLF1);\r
+ end\r
+\r
+ LLF = -1/n * sumLogLLF2;\r
+\r
+ % update distance parameter to check algorithm convergence (delta(phi, Phi))\r
+ deltaPhi = [ deltaPhi, max(max(max((abs(phi-Phi))./(1+abs(phi))))) ];\r
+ if length(deltaPhi) > deltaPhiBufferSize\r
+ deltaPhi = deltaPhi(2:length(deltaPhi));\r
+ end\r
+ sumDeltaPhi = sum(abs(deltaPhi));\r
+\r
+ % update other local variables\r
+ Phi = phi;\r
+ ite = ite+1;\r
+\r
+ end\r
+\r
+end\r