From 23b9fb13bc6e82d7ca43bfb83aa85b6cd69c52c0 Mon Sep 17 00:00:00 2001 From: Benjamin Auder <benjamin.auder@somewhere> Date: Tue, 29 Aug 2017 17:03:37 +0200 Subject: [PATCH] update EMGLLF.c following changes in EMGLLF.R --- pkg/R/EMGLLF.R | 1 - pkg/src/sources/EMGLLF.c | 34 +++++++++++++++++++++------------- 2 files changed, 21 insertions(+), 14 deletions(-) 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); -} +} -- 2.44.0