3 #include <gsl/gsl_linalg.h>
5 // TODO: don't recompute indexes ai(...) and mi(...) when possible
8 const Real
* phiInit
, // parametre initial de moyenne renormalisé
9 const Real
* rhoInit
, // parametre initial de variance renormalisé
10 const Real
* piInit
, // parametre initial des proportions
11 const Real
* gamInit
, // paramètre initial des probabilités a posteriori de chaque échantillon
12 int mini
, // nombre minimal d'itérations dans l'algorithme EM
13 int maxi
, // nombre maximal d'itérations dans l'algorithme EM
14 Real gamma
, // puissance des proportions dans la pénalisation pour un Lasso adaptatif
15 Real lambda
, // valeur du paramètre de régularisation du Lasso
16 const Real
* X
, // régresseurs
17 const Real
* Y
, // réponse
18 Real tau
, // seuil pour accepter la convergence
19 // OUT parameters (all pointers, to be modified)
20 Real
* phi
, // parametre de moyenne renormalisé, calculé par l'EM
21 Real
* rho
, // parametre de variance renormalisé, calculé par l'EM
22 Real
* pi
, // parametre des proportions renormalisé, calculé par l'EM
23 Real
* LLF
, // log vraisemblance associée à cet échantillon, pour les valeurs estimées des paramètres
26 // additional size parameters
27 int n
, // nombre d'echantillons
28 int p
, // nombre de covariables
29 int m
, // taille de Y (multivarié)
30 int k
) // nombre de composantes dans le mélange
33 copyArray(phiInit
, phi
, p
*m
*k
);
34 copyArray(rhoInit
, rho
, m
*m
*k
);
35 copyArray(piInit
, pi
, k
);
37 //S is already allocated, and doesn't need to be 'zeroed'
39 //Other local variables: same as in R
40 Real
* gam
= (Real
*)malloc(n
*k
*sizeof(Real
));
41 copyArray(gamInit
, gam
, n
*k
);
42 Real
* Gram2
= (Real
*)malloc(p
*p
*k
*sizeof(Real
));
43 Real
* ps2
= (Real
*)malloc(p
*m
*k
*sizeof(Real
));
44 Real
* b
= (Real
*)malloc(k
*sizeof(Real
));
45 Real
* X2
= (Real
*)malloc(n
*p
*k
*sizeof(Real
));
46 Real
* Y2
= (Real
*)malloc(n
*m
*k
*sizeof(Real
));
50 Real
* pi2
= (Real
*)malloc(k
*sizeof(Real
));
51 Real
* ps
= (Real
*)malloc(m
*k
*sizeof(Real
));
52 Real
* nY2
= (Real
*)malloc(m
*k
*sizeof(Real
));
53 Real
* ps1
= (Real
*)malloc(n
*m
*k
*sizeof(Real
));
54 Real
* Gam
= (Real
*)malloc(n
*k
*sizeof(Real
));
55 const Real EPS
= 1e-15;
56 // Additional (not at this place, in R file)
57 Real
* gam2
= (Real
*)malloc(k
*sizeof(Real
));
58 Real
* sqNorm2
= (Real
*)malloc(k
*sizeof(Real
));
59 gsl_matrix
* matrix
= gsl_matrix_alloc(m
, m
);
60 gsl_permutation
* permutation
= gsl_permutation_alloc(m
);
61 Real
* YiRhoR
= (Real
*)malloc(m
*sizeof(Real
));
62 Real
* XiPhiR
= (Real
*)malloc(m
*sizeof(Real
));
63 const Real gaussConstM
= pow(2.*M_PI
,m
/2.);
64 Real
* Phi
= (Real
*)malloc(p
*m
*k
*sizeof(Real
));
65 Real
* Rho
= (Real
*)malloc(m
*m
*k
*sizeof(Real
));
66 Real
* Pi
= (Real
*)malloc(k
*sizeof(Real
));
68 while (ite
< mini
|| (ite
< maxi
&& (dist
>= tau
|| dist2
>= sqrt(tau
))))
70 copyArray(phi
, Phi
, p
*m
*k
);
71 copyArray(rho
, Rho
, m
*m
*k
);
74 // Calculs associés a Y et X
75 for (int r
=0; r
<k
; r
++)
77 for (int mm
=0; mm
<m
; mm
++)
79 //Y2[,mm,r] = sqrt(gam[,r]) * Y[,mm]
80 for (int u
=0; u
<n
; u
++)
81 Y2
[ai(u
,mm
,r
,n
,m
,k
)] = sqrt(gam
[mi(u
,r
,n
,k
)]) * Y
[mi(u
,mm
,n
,m
)];
83 for (int i
=0; i
<n
; i
++)
85 //X2[i,,r] = sqrt(gam[i,r]) * X[i,]
86 for (int u
=0; u
<p
; u
++)
87 X2
[ai(i
,u
,r
,n
,p
,k
)] = sqrt(gam
[mi(i
,r
,n
,k
)]) * X
[mi(i
,u
,n
,p
)];
89 for (int mm
=0; mm
<m
; mm
++)
91 //ps2[,mm,r] = crossprod(X2[,,r],Y2[,mm,r])
92 for (int u
=0; u
<p
; u
++)
95 for (int v
=0; v
<n
; v
++)
96 dotProduct
+= X2
[ai(v
,u
,r
,n
,p
,k
)] * Y2
[ai(v
,mm
,r
,n
,m
,k
)];
97 ps2
[ai(u
,mm
,r
,p
,m
,k
)] = dotProduct
;
100 for (int j
=0; j
<p
; j
++)
102 for (int s
=0; s
<p
; s
++)
104 //Gram2[j,s,r] = crossprod(X2[,j,r], X2[,s,r])
105 Real dotProduct
= 0.;
106 for (int u
=0; u
<n
; u
++)
107 dotProduct
+= X2
[ai(u
,j
,r
,n
,p
,k
)] * X2
[ai(u
,s
,r
,n
,p
,k
)];
108 Gram2
[ai(j
,s
,r
,p
,p
,k
)] = dotProduct
;
118 for (int r
=0; r
<k
; r
++)
120 //b[r] = sum(abs(phi[,,r]))
122 for (int u
=0; u
<p
; u
++)
123 for (int v
=0; v
<m
; v
++)
124 sumAbsPhi
+= fabs(phi
[ai(u
,v
,r
,p
,m
,k
)]);
127 //gam2 = colSums(gam)
128 for (int u
=0; u
<k
; u
++)
130 Real sumOnColumn
= 0.;
131 for (int v
=0; v
<n
; v
++)
132 sumOnColumn
+= gam
[mi(v
,u
,n
,k
)];
133 gam2
[u
] = sumOnColumn
;
135 //a = sum(gam %*% log(pi))
137 for (int u
=0; u
<n
; u
++)
139 Real dotProduct
= 0.;
140 for (int v
=0; v
<k
; v
++)
141 dotProduct
+= gam
[mi(u
,v
,n
,k
)] * log(pi
[v
]);
145 //tant que les proportions sont negatives
147 int pi2AllPositive
= 0;
149 while (!pi2AllPositive
)
151 //pi2 = pi + 0.1^kk * ((1/n)*gam2 - pi)
152 Real pow_01_kk
= pow(0.1,kk
);
153 for (int r
=0; r
<k
; r
++)
154 pi2
[r
] = pi
[r
] + pow_01_kk
* (invN
*gam2
[r
] - pi
[r
]);
155 //pi2AllPositive = all(pi2 >= 0)
157 for (int r
=0; r
<k
; r
++)
169 Real piPowGammaDotB
= 0.;
170 for (int v
=0; v
<k
; v
++)
171 piPowGammaDotB
+= pow(pi
[v
],gamma
) * b
[v
];
173 Real pi2PowGammaDotB
= 0.;
174 for (int v
=0; v
<k
; v
++)
175 pi2PowGammaDotB
+= pow(pi2
[v
],gamma
) * b
[v
];
176 //sum(gam2 * log(pi2))
177 Real gam2DotLogPi2
= 0.;
178 for (int v
=0; v
<k
; v
++)
179 gam2DotLogPi2
+= gam2
[v
] * log(pi2
[v
]);
181 //t(m) la plus grande valeur dans la grille O.1^k tel que ce soit décroissante ou constante
182 while (-invN
*a
+ lambda
*piPowGammaDotB
< -invN
*gam2DotLogPi2
+ lambda
*pi2PowGammaDotB
185 Real pow_01_kk
= pow(0.1,kk
);
186 //pi2 = pi + 0.1^kk * (1/n*gam2 - pi)
187 for (int v
=0; v
<k
; v
++)
188 pi2
[v
] = pi
[v
] + pow_01_kk
* (invN
*gam2
[v
] - pi
[v
]);
189 //pi2 was updated, so we recompute pi2PowGammaDotB and gam2DotLogPi2
190 pi2PowGammaDotB
= 0.;
191 for (int v
=0; v
<k
; v
++)
192 pi2PowGammaDotB
+= pow(pi2
[v
],gamma
) * b
[v
];
194 for (int v
=0; v
<k
; v
++)
195 gam2DotLogPi2
+= gam2
[v
] * log(pi2
[v
]);
198 Real t
= pow(0.1,kk
);
199 //sum(pi + t*(pi2-pi))
200 Real sumPiPlusTbyDiff
= 0.;
201 for (int v
=0; v
<k
; v
++)
202 sumPiPlusTbyDiff
+= (pi
[v
] + t
*(pi2
[v
] - pi
[v
]));
203 //pi = (pi + t*(pi2-pi)) / sum(pi + t*(pi2-pi))
204 for (int v
=0; v
<k
; v
++)
205 pi
[v
] = (pi
[v
] + t
*(pi2
[v
] - pi
[v
])) / sumPiPlusTbyDiff
;
208 for (int r
=0; r
<k
; r
++)
210 for (int mm
=0; mm
<m
; mm
++)
212 for (int i
=0; i
<n
; i
++)
214 //< X2[i,,r] , phi[,mm,r] >
215 Real dotProduct
= 0.;
216 for (int u
=0; u
<p
; u
++)
217 dotProduct
+= X2
[ai(i
,u
,r
,n
,p
,k
)] * phi
[ai(u
,mm
,r
,p
,m
,k
)];
218 //ps1[i,mm,r] = Y2[i,mm,r] * sum(X2[i,,r] * phi[,mm,r])
219 ps1
[ai(i
,mm
,r
,n
,m
,k
)] = Y2
[ai(i
,mm
,r
,n
,m
,k
)] * dotProduct
;
221 //ps[mm,r] = sum(ps1[,mm,r])
223 for (int u
=0; u
<n
; u
++)
224 sumPs1
+= ps1
[ai(u
,mm
,r
,n
,m
,k
)];
225 ps
[mi(mm
,r
,m
,k
)] = sumPs1
;
226 //nY2[mm,r] = sum(Y2[,mm,r]^2)
228 for (int u
=0; u
<n
; u
++)
229 sumY2
+= Y2
[ai(u
,mm
,r
,n
,m
,k
)] * Y2
[ai(u
,mm
,r
,n
,m
,k
)];
230 nY2
[mi(mm
,r
,m
,k
)] = sumY2
;
231 //rho[mm,mm,r] = (ps[mm,r]+sqrt(ps[mm,r]^2+4*nY2[mm,r]*(gam2[r]))) / (2*nY2[mm,r])
232 rho
[ai(mm
,mm
,r
,m
,m
,k
)] = ( ps
[mi(mm
,r
,m
,k
)] + sqrt( ps
[mi(mm
,r
,m
,k
)]*ps
[mi(mm
,r
,m
,k
)]
233 + 4*nY2
[mi(mm
,r
,m
,k
)] * gam2
[r
] ) ) / (2*nY2
[mi(mm
,r
,m
,k
)]);
237 for (int r
=0; r
<k
; r
++)
239 for (int j
=0; j
<p
; j
++)
241 for (int mm
=0; mm
<m
; mm
++)
243 //sum(phi[-j,mm,r] * Gram2[j, setdiff(1:p,j),r])
244 Real phiDotGram2
= 0.;
245 for (int u
=0; u
<p
; u
++)
248 phiDotGram2
+= phi
[ai(u
,mm
,r
,p
,m
,k
)] * Gram2
[ai(j
,u
,r
,p
,p
,k
)];
250 //S[j,mm,r] = -rho[mm,mm,r]*ps2[j,mm,r] + sum(phi[-j,mm,r] * Gram2[j,-j,r])
251 S
[ai(j
,mm
,r
,p
,m
,k
)] = -rho
[ai(mm
,mm
,r
,m
,m
,k
)] * ps2
[ai(j
,mm
,r
,p
,m
,k
)] + phiDotGram2
;
252 Real pirPowGamma
= pow(pi
[r
],gamma
);
253 if (fabs(S
[ai(j
,mm
,r
,p
,m
,k
)]) <= n
*lambda
*pirPowGamma
)
254 phi
[ai(j
,mm
,r
,p
,m
,k
)] = 0.;
255 else if (S
[ai(j
,mm
,r
,p
,m
,k
)] > n
*lambda
*pirPowGamma
)
257 phi
[ai(j
,mm
,r
,p
,m
,k
)] = (n
*lambda
*pirPowGamma
- S
[ai(j
,mm
,r
,p
,m
,k
)])
258 / Gram2
[ai(j
,j
,r
,p
,p
,k
)];
262 phi
[ai(j
,mm
,r
,p
,m
,k
)] = -(n
*lambda
*pirPowGamma
+ S
[ai(j
,mm
,r
,p
,m
,k
)])
263 / Gram2
[ai(j
,j
,r
,p
,p
,k
)];
274 Real sumLogLLF2
= 0.;
275 for (int i
=0; i
<n
; i
++)
277 for (int r
=0; r
<k
; r
++)
279 //compute Y[i,]%*%rho[,,r]
280 for (int u
=0; u
<m
; u
++)
283 for (int v
=0; v
<m
; v
++)
284 YiRhoR
[u
] += Y
[mi(i
,v
,n
,m
)] * rho
[ai(v
,u
,r
,m
,m
,k
)];
287 //compute X[i,]%*%phi[,,r]
288 for (int u
=0; u
<m
; u
++)
291 for (int v
=0; v
<p
; v
++)
292 XiPhiR
[u
] += X
[mi(i
,v
,n
,p
)] * phi
[ai(v
,u
,r
,p
,m
,k
)];
295 //compute sq norm || Y(:,i)*rho(:,:,r)-X(i,:)*phi(:,:,r) ||_2^2
297 for (int u
=0; u
<m
; u
++)
298 sqNorm2
[r
] += (YiRhoR
[u
]-XiPhiR
[u
]) * (YiRhoR
[u
]-XiPhiR
[u
]);
303 for (int r
=0; r
<k
; r
++)
305 //compute det(rho[,,r]) [TODO: avoid re-computations]
306 for (int u
=0; u
<m
; u
++)
308 for (int v
=0; v
<m
; v
++)
309 matrix
->data
[u
*m
+v
] = rho
[ai(u
,v
,r
,m
,m
,k
)];
311 gsl_linalg_LU_decomp(matrix
, permutation
, &signum
);
312 Real detRhoR
= gsl_linalg_LU_det(matrix
, signum
);
313 Gam
[mi(i
,r
,n
,k
)] = pi
[r
] * exp(-.5*sqNorm2
[r
]) * detRhoR
;
314 sumLLF1
+= Gam
[mi(i
,r
,n
,k
)] / gaussConstM
;
315 sumGamI
+= Gam
[mi(i
,r
,n
,k
)];
318 sumLogLLF2
+= log(sumLLF1
);
319 for (int r
=0; r
<k
; r
++)
321 //gam[i,] = Gam[i,] / sumGamI
322 gam
[mi(i
,r
,n
,k
)] = sumGamI
> EPS
? Gam
[mi(i
,r
,n
,k
)] / sumGamI
: 0.;
326 //sumPen = sum(pi^gamma * b)
328 for (int r
=0; r
<k
; r
++)
329 sumPen
+= pow(pi
[r
],gamma
) * b
[r
];
330 //LLF[ite] = -sumLogLLF2/n + lambda*sumPen
331 LLF
[ite
] = -invN
* sumLogLLF2
+ lambda
* sumPen
;
332 dist
= ite
==0 ? LLF
[ite
] : (LLF
[ite
] - LLF
[ite
-1]) / (1. + fabs(LLF
[ite
]));
334 //Dist1 = max( abs(phi-Phi) / (1+abs(phi)) )
336 for (int u
=0; u
<p
; u
++)
338 for (int v
=0; v
<m
; v
++)
340 for (int w
=0; w
<k
; w
++)
342 Real tmpDist
= fabs(phi
[ai(u
,v
,w
,p
,m
,k
)]-Phi
[ai(u
,v
,w
,p
,m
,k
)])
343 / (1.+fabs(phi
[ai(u
,v
,w
,p
,m
,k
)]));
349 //Dist2 = max( (abs(rho-Rho)) / (1+abs(rho)) )
351 for (int u
=0; u
<m
; u
++)
353 for (int v
=0; v
<m
; v
++)
355 for (int w
=0; w
<k
; w
++)
357 Real tmpDist
= fabs(rho
[ai(u
,v
,w
,m
,m
,k
)]-Rho
[ai(u
,v
,w
,m
,m
,k
)])
358 / (1.+fabs(rho
[ai(u
,v
,w
,m
,m
,k
)]));
364 //Dist3 = max( (abs(pi-Pi)) / (1+abs(Pi)))
366 for (int u
=0; u
<n
; u
++)
368 for (int v
=0; v
<k
; v
++)
370 Real tmpDist
= fabs(pi
[v
]-Pi
[v
]) / (1.+fabs(pi
[v
]));
375 //dist2=max([max(Dist1),max(Dist2),max(Dist3)]);
385 //affec = apply(gam, 1, which.max)
386 for (int i
=0; i
<n
; i
++)
390 for (int j
=0; j
<k
; j
++)
392 if (gam
[mi(i
,j
,n
,k
)] > rowMax
)
394 affec
[i
] = j
+1; //R indices start at 1
395 rowMax
= gam
[mi(i
,j
,n
,k
)];
412 gsl_matrix_free(matrix
);
413 gsl_permutation_free(permutation
);