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