prepare structure for R package

[valse.git] / OLD_MATLAB / ProcLassoRank / EMGrank.m

diff --git a/OLD_MATLAB/ProcLassoRank/EMGrank.m b/OLD_MATLAB/ProcLassoRank/EMGrank.m
new file mode 100644 (file)
index 0000000..6ffc313
--- /dev/null
+++ b/OLD_MATLAB/ProcLassoRank/EMGrank.m
@@ -0,0 +1,69 @@
+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