Inputparametters.R + SelecModel
[valse.git] / ProcLassoMLE / EMGLLF.m
1 function[phi,rho,pi,LLF,S] = EMGLLF(phiInit,rhoInit,piInit,gamInit,mini,maxi,gamma,lambda,X,Y,tau)
2
3 %Get matrices dimensions
4 PI = 4.0 * atan(1.0);
5 n = size(X, 1);
6 [p,m,k] = size(phiInit);
7
8 %Initialize outputs
9 phi = phiInit;
10 rho = rhoInit;
11 pi = piInit;
12 LLF = zeros(maxi,1);
13 S = zeros(p,m,k);
14
15 %Other local variables
16 %NOTE: variables order is always n,p,m,k
17 gam = gamInit;
18 Gram2 = zeros(p,p,k);
19 ps2 = zeros(p,m,k);
20 b = zeros(k,1);
21 pen = zeros(maxi,k);
22 X2 = zeros(n,p,k);
23 Y2 = zeros(n,m,k);
24 dist = 0;
25 dist2 = 0;
26 ite = 1;
27 pi2 = zeros(k,1);
28 ps = zeros(m,k);
29 nY2 = zeros(m,k);
30 ps1 = zeros(n,m,k);
31 nY21 = zeros(n,m,k);
32 Gam = zeros(n,k);
33 EPS = 1e-15;
34
35 while ite<=mini || (ite<=maxi && (dist>=tau || dist2>=sqrt(tau)))
36
37 Phi = phi;
38 Rho = rho;
39 Pi = pi;
40
41 %Calculs associés à Y et X
42 for r=1:k
43 for mm=1:m
44 Y2(:,mm,r) = sqrt(gam(:,r)) .* Y(:,mm);
45 end
46 for i=1:n
47 X2(i,:,r) = X(i,:) .* sqrt(gam(i,r));
48 end
49 for mm=1:m
50 ps2(:,mm,r) = transpose(X2(:,:,r)) * Y2(:,mm,r);
51 end
52 for j=1:p
53 for s=1:p
54 Gram2(j,s,r) = dot(X2(:,j,r), X2(:,s,r));
55 end
56 end
57 end
58
59 %%%%%%%%%%
60 %Etape M %
61 %%%%%%%%%%
62
63 %Pour pi
64 for r=1:k
65 b(r) = sum(sum(abs(phi(:,:,r))));
66 end
67 gam2 = sum(gam,1);
68 a = sum(gam*transpose(log(pi)));
69
70 %tant que les proportions sont negatives
71 kk = 0;
72 pi2AllPositive = false;
73 while ~pi2AllPositive
74 pi2 = pi + 0.1^kk * ((1/n)*gam2 - pi);
75 pi2AllPositive = true;
76 for r=1:k
77 if pi2(r) < 0
78 pi2AllPositive = false;
79 break;
80 end
81 end
82 kk = kk+1;
83 end
84
85 %t(m) la plus grande valeur dans la grille O.1^k tel que ce soit
86 %décroissante ou constante
87 while (-1/n*a+lambda*((pi.^gamma)*b))<(-1/n*gam2*transpose(log(pi2))+lambda.*(pi2.^gamma)*b) && kk<1000
88 pi2 = pi+0.1^kk*(1/n*gam2-pi);
89 kk = kk+1;
90 end
91 t = 0.1^(kk);
92 pi = (pi+t*(pi2-pi)) / sum(pi+t*(pi2-pi));
93
94 %Pour phi et rho
95 for r=1:k
96 for mm=1:m
97 for i=1:n
98 ps1(i,mm,r) = Y2(i,mm,r) * dot(X2(i,:,r), phi(:,mm,r));
99 nY21(i,mm,r) = (Y2(i,mm,r))^2;
100 end
101 ps(mm,r) = sum(ps1(:,mm,r));
102 nY2(mm,r) = sum(nY21(:,mm,r));
103 rho(mm,mm,r) = ((ps(mm,r)+sqrt(ps(mm,r)^2+4*nY2(mm,r)*(gam2(r))))/(2*nY2(mm,r)));
104 end
105 end
106 for r=1:k
107 for j=1:p
108 for mm=1:m
109 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)')...
110 + dot(phi(j+1:p,mm,r),Gram2(j,j+1:p,r)');
111 if abs(S(j,mm,r)) <= n*lambda*(pi(r)^gamma)
112 phi(j,mm,r)=0;
113 else
114 if S(j,mm,r)> n*lambda*(pi(r)^gamma)
115 phi(j,mm,r)=(n*lambda*(pi(r)^gamma)-S(j,mm,r))/Gram2(j,j,r);
116 else
117 phi(j,mm,r)=-(n*lambda*(pi(r)^gamma)+S(j,mm,r))/Gram2(j,j,r);
118 end
119 end
120 end
121 end
122 end
123
124 %%%%%%%%%%
125 %Etape E %
126 %%%%%%%%%%
127
128 sumLogLLF2 = 0.0;
129 for i=1:n
130 %precompute dot products to numerically adjust their values
131 dotProducts = zeros(k,1);
132 for r=1:k
133 dotProducts(r)= (Y(i,:)*rho(:,:,r)-X(i,:)*phi(:,:,r)) * transpose(Y(i,:)*rho(:,:,r)-X(i,:)*phi(:,:,r));
134 end
135 shift = 0.5*min(dotProducts);
136
137 %compute Gam(:,:) using shift determined above
138 sumLLF1 = 0.0;
139 for r=1:k
140 Gam(i,r) = pi(r)*det(rho(:,:,r))*exp(-0.5*dotProducts(r) + shift);
141 sumLLF1 = sumLLF1 + Gam(i,r)/(2*PI)^(m/2);
142 end
143 sumLogLLF2 = sumLogLLF2 + log(sumLLF1);
144 sumGamI = sum(Gam(i,:));
145 if sumGamI > EPS
146 gam(i,:) = Gam(i,:) / sumGamI;
147 else
148 gam(i,:) = zeros(k,1);
149 end
150 end
151
152 sumPen = 0.0;
153 for r=1:k
154 sumPen = sumPen + pi(r).^gamma .* b(r);
155 end
156 LLF(ite) = -(1/n)*sumLogLLF2 + lambda*sumPen;
157
158 if ite == 1
159 dist = LLF(ite);
160 else
161 dist = (LLF(ite)-LLF(ite-1))/(1+abs(LLF(ite)));
162 end
163
164 Dist1=max(max(max((abs(phi-Phi))./(1+abs(phi)))));
165 Dist2=max(max(max((abs(rho-Rho))./(1+abs(rho)))));
166 Dist3=max(max((abs(pi-Pi))./(1+abs(Pi))));
167 dist2=max([Dist1,Dist2,Dist3]);
168
169 ite=ite+1;
170 end
171
172 pi = transpose(pi);
173
174 end