#include <gsl/gsl_linalg.h>
// TODO: don't recompute indexes ai(...) and mi(...) when possible
-void EMGLLH_core(
+void EMGLLF_core(
// IN parameters
const Real* phiInit, // parametre initial de moyenne renormalisé
const Real* rhoInit, // parametre initial de variance renormalisé
Real lambda, // valeur du paramètre de régularisation du Lasso
const Real* X, // régresseurs
const Real* Y, // réponse
- Real tau, // seuil pour accepter la convergence
+ Real eps, // seuil pour accepter la convergence
// OUT parameters (all pointers, to be modified)
Real* phi, // parametre de moyenne renormalisé, calculé par l'EM
Real* rho, // parametre de variance renormalisé, calculé par l'EM
//Other local variables: same as in R
Real* gam = (Real*)malloc(n*k*sizeof(Real));
+ Real* logGam = (Real*)malloc(k*sizeof(Real));
copyArray(gamInit, gam, n*k);
Real* Gram2 = (Real*)malloc(p*p*k*sizeof(Real));
Real* ps2 = (Real*)malloc(p*m*k*sizeof(Real));
Real* Y2 = (Real*)malloc(n*m*k*sizeof(Real));
*llh = -INFINITY;
Real* pi2 = (Real*)malloc(k*sizeof(Real));
- const Real EPS = 1e-15;
// Additional (not at this place, in R file)
Real* gam2 = (Real*)malloc(k*sizeof(Real));
Real* sqNorm2 = (Real*)malloc(k*sizeof(Real));
Real* Rho = (Real*)malloc(m*m*k*sizeof(Real));
Real* Pi = (Real*)malloc(k*sizeof(Real));
- for (int ite=0; ite<maxi; ite++)
+ for (int ite=1; ite<=maxi; ite++)
{
copyArray(phi, Phi, p*m*k);
copyArray(rho, Rho, m*m*k);
/////////////
// Precompute det(rho[,,r]) for r in 1...k
+ int signum;
for (int r=0; r<k; r++)
{
for (int u=0; u<m; u++)
detRho[r] = gsl_linalg_LU_det(matrix, signum);
}
- int signum;
Real sumLogLLH = 0.;
for (int i=0; i<n; i++)
{
sqNorm2[r] += (YiRhoR[u]-XiPhiR[u]) * (YiRhoR[u]-XiPhiR[u]);
}
- Real sumGamI = 0.;
+ // Update gam[,]; use log to avoid numerical problems
+ Real maxLogGam = -INFINITY;
for (int r=0; r<k; r++)
{
- gam[mi(i,r,n,k)] = pi[r] * exp(-.5*sqNorm2[r]) * detRho[r];
- sumGamI += gam[mi(i,r,n,k)];
+ logGam[r] = log(pi[r]) - .5 * sqNorm2[r] + log(detRho[r]);
+ if (maxLogGam < logGam[r])
+ maxLogGam = logGam[r];
}
-
- sumLogLLH += log(sumGamI) - log(gaussConstM);
- if (sumGamI > EPS) //else: gam[i,] is already ~=0
+ Real norm_fact = 0.;
+ for (int r=0; r<k; r++)
{
- for (int r=0; r<k; r++)
- gam[mi(i,r,n,k)] /= sumGamI;
+ logGam[r] = logGam[r] - maxLogGam; //adjust without changing proportions
+ gam[mi(i,r,n,k)] = exp(logGam[r]); //gam[i, ] <- exp(logGam)
+ norm_fact += gam[mi(i,r,n,k)]; //norm_fact <- sum(gam[i, ])
}
+ // gam[i, ] <- gam[i, ] / norm_fact
+ for (int r=0; r<k; r++)
+ gam[mi(i,r,n,k)] /= norm_fact;
+
+ sumLogLLH += log(norm_fact) - log(gaussConstM);
}
//sumPen = sum(pi^gamma * b)
for (int r=0; r<k; r++)
sumPen += pow(pi[r],gamma) * b[r];
Real last_llh = *llh;
- //llh = -sumLogLLH/n + lambda*sumPen
- *llh = -invN * sumLogLLH + lambda * sumPen;
- Real dist = ite==0 ? *llh : (*llh - last_llh) / (1. + fabs(*llh));
+ //llh = -sumLogLLH/n #+ lambda*sumPen
+ *llh = -invN * sumLogLLH; //+ lambda * sumPen;
+ Real dist = ( ite==1 ? *llh : (*llh - last_llh) / (1. + fabs(*llh)) );
//Dist1 = max( abs(phi-Phi) / (1+abs(phi)) )
Real Dist1 = 0.;
if (Dist3 > dist2)
dist2 = Dist3;
- if (ite >= mini && (dist >= tau || dist2 >= sqrt(tau)))
+ if (ite >= mini && (dist >= eps || dist2 >= sqrt(eps)))
break;
}
//free memory
free(b);
free(gam);
+ free(logGam);
free(Phi);
free(Rho);
free(Pi);
free(Gram2);
free(ps2);
+ free(detRho);
gsl_matrix_free(matrix);
gsl_permutation_free(permutation);
free(XiPhiR);