reorganize test folder: tests generation with MATLAB almost OK (but suboptimal method)
[valse.git] / src / test / generate_test_data / MATLAB_test_helpers / EMGrank.m
... / ...
CommitLineData
1function[phi,LLF] = EMGrank(Pi,Rho,mini,maxi,X,Y,tau,rank)\r
2\r
3 % get matrix sizes\r
4 [~,m,k] = size(Rho);\r
5 [n,p] = size(X);\r
6\r
7 % allocate output matrices\r
8 phi = zeros(p,m,k);\r
9 Z = ones(n,1,'int64');\r
10 LLF = 0.0;\r
11\r
12 % local variables\r
13 Phi = zeros(p,m,k);\r
14 deltaPhi = 0.0;\r
15 deltaPhi = [];\r
16 sumDeltaPhi = 0.0;\r
17 deltaPhiBufferSize = 20;\r
18\r
19 %main loop (at least mini iterations)\r
20 ite = int64(1);\r
21 while ite<=mini || (ite<=maxi && sumDeltaPhi>tau)\r
22\r
23 %M step: Mise à jour de Beta (et donc phi)\r
24 for r=1:k\r
25 if (sum(Z==r) == 0)\r
26 continue;\r
27 end\r
28 %U,S,V = SVD of (t(Xr)Xr)^{-1} * t(Xr) * Yr\r
29 [U,S,V] = svd(pinv(transpose(X(Z==r,:))*X(Z==r,:))*transpose(X(Z==r,:))*Y(Z==r,:));\r
30 %Set m-rank(r) singular values to zero, and recompose\r
31 %best rank(r) approximation of the initial product\r
32 S(rank(r)+1:end,:) = 0;\r
33 phi(:,:,r) = U * S * transpose(V) * Rho(:,:,r);\r
34 end\r
35\r
36 %Etape E et calcul de LLF\r
37 sumLogLLF2 = 0.0;\r
38 for i=1:n\r
39 sumLLF1 = 0.0;\r
40 maxLogGamIR = -Inf;\r
41 for r=1:k\r
42 dotProduct = (Y(i,:)*Rho(:,:,r)-X(i,:)*phi(:,:,r)) * transpose(Y(i,:)*Rho(:,:,r)-X(i,:)*phi(:,:,r));\r
43 logGamIR = log(Pi(r)) + log(det(Rho(:,:,r))) - 0.5*dotProduct;\r
44 %Z(i) = index of max (gam(i,:))\r
45 if logGamIR > maxLogGamIR\r
46 Z(i) = r;\r
47 maxLogGamIR = logGamIR;\r
48 end\r
49 sumLLF1 = sumLLF1 + exp(logGamIR) / (2*pi)^(m/2);\r
50 end\r
51 sumLogLLF2 = sumLogLLF2 + log(sumLLF1);\r
52 end\r
53\r
54 LLF = -1/n * sumLogLLF2;\r
55\r
56 % update distance parameter to check algorithm convergence (delta(phi, Phi))\r
57 deltaPhi = [ deltaPhi, max(max(max((abs(phi-Phi))./(1+abs(phi))))) ];\r
58 if length(deltaPhi) > deltaPhiBufferSize\r
59 deltaPhi = deltaPhi(2:length(deltaPhi));\r
60 end\r
61 sumDeltaPhi = sum(abs(deltaPhi));\r
62\r
63 % update other local variables\r
64 Phi = phi;\r
65 ite = ite+1;\r
66\r
67 end\r
68\r
69end\r