From: Benjamin Auder <benjamin.auder@somewhere>
Date: Tue, 29 Aug 2017 15:03:37 +0000 (+0200)
Subject: update EMGLLF.c following changes in EMGLLF.R
X-Git-Url: https://git.auder.net/%7B%7B%20asset%28%27mixstore/images/scripts/doc/pieces/%7B%7B?a=commitdiff_plain;h=23b9fb13bc6e82d7ca43bfb83aa85b6cd69c52c0;p=valse.git

update EMGLLF.c following changes in EMGLLF.R
---

diff --git a/pkg/R/EMGLLF.R b/pkg/R/EMGLLF.R
index 2aeea53..0d8607c 100644
--- a/pkg/R/EMGLLF.R
+++ b/pkg/R/EMGLLF.R
@@ -74,7 +74,6 @@ EMGLLF <- function(phiInit, rhoInit, piInit, gamInit, mini, maxi, gamma, lambda,
   ps2 <- array(0, dim = c(p, m, k))
   X2 <- array(0, dim = c(n, p, k))
   Y2 <- array(0, dim = c(n, m, k))
-  EPS <- 1e-15
 
   for (ite in 1:maxi)
   {
diff --git a/pkg/src/sources/EMGLLF.c b/pkg/src/sources/EMGLLF.c
index d2f5a8e..d9f92c0 100644
--- a/pkg/src/sources/EMGLLF.c
+++ b/pkg/src/sources/EMGLLF.c
@@ -39,6 +39,7 @@ void EMGLLF_core(
 
 	//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));
@@ -47,7 +48,6 @@ void EMGLLF_core(
 	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));
@@ -300,19 +300,26 @@ void EMGLLF_core(
 					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)
@@ -320,9 +327,9 @@ void EMGLLF_core(
 		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==1 ? *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.;
@@ -394,6 +401,7 @@ void EMGLLF_core(
 	//free memory
 	free(b);
 	free(gam);
+	free(logGam);
 	free(Phi);
 	free(Rho);
 	free(Pi);
@@ -409,4 +417,4 @@ void EMGLLF_core(
 	free(X2);
 	free(Y2);
 	free(sqNorm2);
-}
+}