Changeset 3648 for trunk/test/poisson.cc

Timestamp:

Jun 2, 2017, 8:01:27 AM (6 years ago)

Author:

Peter

Message:

add some tests for regression::Poisson. refs #882

File:

• trunk/test/poisson.cc (modified) (3 diffs)

trunk/test/poisson.cc

@@ -2 +2 @@
 
 /*
-  Copyright (C) 2017 Jari Häkkinen
+  Copyright (C) 2017 Peter Johansson
 
   This file is part of the yat library, http://dev.thep.lu.se/yat
@@ -25 +25 @@
 
 #include "yat/regression/Poisson.h"
+#include "yat/random/random.h"
+#include "yat/statistics/Averager.h"
 #include "yat/utility/Matrix.h"
 #include "yat/utility/Vector.h"
 
+#include <vector>
+
 using namespace theplu::yat;
+
+void analyse(const utility::Vector& b,
+             std::vector<statistics::Averager>& stats,
+             std::vector<statistics::Averager>& stats2);
 
 int main(int argc, char* argv[])
@@ -48 +56 @@
   model.predict(X.row_const_view(0));
 
+  utility::Vector b(4);
+  b(0) = 0.5;
+  b(1) = 2.0;
+  b(2) = 0.75;
+  b(3) = -1.25;
+  std::vector<statistics::Averager> stats(b.size());
+  std::vector<statistics::Averager> stats2(b.size());
+  for (size_t i=0; i<100; ++i)
+    analyse(b, stats, stats2);
+
+  for (size_t i=0; i<b.size(); ++i) {
+    suite.out() << i << " " << b(i) << " " << stats[i].mean() << " "
+                << stats[i].standard_error() << "\n";
+    if (stats[i].standard_error() == 0.0) {
+      suite.xadd(false);
+      suite.err() << "error: standard error is 0.0\n";
+    }
+    else if (std::abs(stats[i].mean() - b(i)) > 5*stats[i].standard_error()) {
+      suite.xadd(false);
+      suite.err() << "error: average for param " << i << ": "
+                  << stats[i].mean() << "\n";
+    }
+  }
+
   return suite.return_value();
 }
+
+
+void generate_data(const utility::Vector& b,
+                   utility::Matrix& X, utility::Vector& y)
+{
+  size_t n = 5000;
+  X.resize(n, b.size()-1);
+  y.resize(n);
+  random::Gaussian gauss;
+  random::Poisson poisson;
+  for (utility::Matrix::iterator it=X.begin(); it!=X.end(); ++it)
+    *it = gauss();
+
+  for (size_t i=0; i<n; ++i) {
+    double lnmu = b(0);
+    for (size_t j=1; j<b.size(); ++j)
+      lnmu += X(i, j-1) * b(j);
+
+    double mu = exp(lnmu);
+    y(i) = poisson(mu);
+  }
+}
+
+
+void analyse(const utility::Vector& b,
+             std::vector<statistics::Averager>& stats,
+             std::vector<statistics::Averager>& stats2)
+{
+  utility::Matrix X;
+  utility::Vector y;
+  generate_data(b, X, y);
+  theplu::yat::regression::Poisson poisson;
+  poisson.fit(X, y);
+  for (size_t i=0; i<poisson.fit_parameters().size(); ++i) {
+    stats[i].add(poisson.fit_parameters()(i));
+  }
+}