[Genabel-commits] r1661 - pkg/ProbABEL/src

[noreply at r-forge.r-project.org]

Author: lckarssen
Date: 2014-03-27 10:40:53 +0100 (Thu, 27 Mar 2014)
New Revision: 1661

Modified:
   pkg/ProbABEL/src/reg1.cpp
Log:
Reverting the log-of-products part of commit r1658, because it can lead to underflows, as Yurii pointed out.


Modified: pkg/ProbABEL/src/reg1.cpp
===================================================================
--- pkg/ProbABEL/src/reg1.cpp	2014-03-25 21:24:03 UTC (rev 1660)
+++ pkg/ProbABEL/src/reg1.cpp	2014-03-27 09:40:53 UTC (rev 1661)
@@ -780,23 +780,12 @@
         }
         // std::cout << "beta:\n"; beta.print();
 
-        // Compute the likelihood. The following commented code gives
-        // the 'easy to understand' algorithm. The code that's
-        // actually used is mathematically equivalent (remember:
-        // log(a*b) = log(a)+log(b)), but faster because log() is
-        // relatively expensive.
-        //    for (int i = 0; i < eMu.nrow; i++) {
-        //          loglik += rdata.Y[i] * eMu_us[i] - log(1. +
-        //                    exp(eMu_us[i]));
-        //    }
+        // Compute the likelihood.
         prevlik = loglik;
-        loglik = 0.;
-        double logterm = 1;
+        loglik = 0;
         for (int i = 0; i < eMu.nrow; i++) {
-            loglik  += rdata.Y[i] * eMu_us[i];
-            logterm *= 1. + exp(eMu_us[i]);
+            loglik += rdata.Y[i] * eMu_us[i] - log(1. + exp(eMu_us[i]));
         }
-        loglik += - log(logterm);
 
         delta = fabs(1. - (prevlik / loglik));
         niter++;