[Genabel-commits] r1531 - branches/ProbABEL-0.50/src

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

Author: maartenk
Date: 2014-01-10 21:21:52 +0100 (Fri, 10 Jan 2014)
New Revision: 1531

Modified:
   branches/ProbABEL-0.50/src/reg1.cpp
Log:
added calculation of loglikl. in EIGEN specific code

Modified: branches/ProbABEL-0.50/src/reg1.cpp
===================================================================
--- branches/ProbABEL-0.50/src/reg1.cpp	2014-01-10 16:02:27 UTC (rev 1530)
+++ branches/ProbABEL-0.50/src/reg1.cpp	2014-01-10 20:21:52 UTC (rev 1531)
@@ -485,14 +485,26 @@
     //cout << endl;
     loglik = 0.;
     double halfrecsig2 = .5 / sigma2;
+#if EIGEN
+    double intercept = beta.get(0, 0);
+    residuals.data= rdata.Y.data.array()-intercept;
+    Eigen::ArrayXXd betacol = beta.data.block(1,0,beta.data.rows()-1,1).array().transpose();
+    Eigen::ArrayXXd resid_sub = (X.data.block(0,1,X.data.rows(),X.data.cols()-1)*betacol.matrix().asDiagonal()).rowwise().sum() ;
+    residuals.data-=resid_sub.matrix();
+    loglik-=(residuals.data.array().square()*halfrecsig2).sum();
+
+#else
     for (int i = 0; i < rdata.nids; i++)
-    {
-        double resid = rdata.Y[i] - beta.get(0, 0); // intercept
-        for (int j = 1; j < beta.nrow; j++)
-            resid -= beta.get(j, 0) * X.get(i, j);
-        residuals[i] = resid;
-        loglik -= halfrecsig2 * resid * resid;
-    }
+     {
+         double resid = rdata.Y[i] - beta.get(0, 0); // intercept
+         for (int j = 1; j < beta.nrow; j++){
+             resid -= beta.get(j, 0) * X.get(i, j);
+
+         }
+         residuals[i] = resid;
+         loglik -= halfrecsig2 * resid * resid;
+     }
+#endif
     loglik -= static_cast<double>(rdata.nids) * log(sqrt(sigma2));
     // cout << "estimate " << rdata.nids << "\n";
     //