source: trunk/test/poisson.cc @ 3648

Last change on this file since 3648 was 3648, checked in by Peter, 6 years ago

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

  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 3.1 KB
Line 
1// $Id: poisson.cc 3648 2017-06-02 06:01:27Z peter $
2
3/*
4  Copyright (C) 2017 Peter Johansson
5
6  This file is part of the yat library, http://dev.thep.lu.se/yat
7
8  The yat library is free software; you can redistribute it and/or
9  modify it under the terms of the GNU General Public License as
10  published by the Free Software Foundation; either version 3 of the
11  License, or (at your option) any later version.
12
13  The yat library is distributed in the hope that it will be useful,
14  but WITHOUT ANY WARRANTY; without even the implied warranty of
15  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
16  General Public License for more details.
17
18  You should have received a copy of the GNU General Public License
19  along with yat. If not, see <http://www.gnu.org/licenses/>.
20*/
21
22#include <config.h>
23
24#include "Suite.h"
25
26#include "yat/regression/Poisson.h"
27#include "yat/random/random.h"
28#include "yat/statistics/Averager.h"
29#include "yat/utility/Matrix.h"
30#include "yat/utility/Vector.h"
31
32#include <vector>
33
34using namespace theplu::yat;
35
36void analyse(const utility::Vector& b,
37             std::vector<statistics::Averager>& stats,
38             std::vector<statistics::Averager>& stats2);
39
40int main(int argc, char* argv[])
41{
42  test::Suite suite(argc, argv);
43
44  regression::Poisson model;
45  size_t n = 20;
46  size_t p = 4;
47  utility::Vector y(n);
48  utility::Matrix X(n, p);
49  model.fit(X, y);
50  if (model.fit_parameters().size() != p+1) {
51    suite.add(false);
52    suite.err() << "error: size of fit_parameters: "
53                << model.fit_parameters().size()
54                << " expected " << p+1 << "\n";
55  }
56  model.predict(X.row_const_view(0));
57
58  utility::Vector b(4);
59  b(0) = 0.5;
60  b(1) = 2.0;
61  b(2) = 0.75;
62  b(3) = -1.25;
63  std::vector<statistics::Averager> stats(b.size());
64  std::vector<statistics::Averager> stats2(b.size());
65  for (size_t i=0; i<100; ++i)
66    analyse(b, stats, stats2);
67
68  for (size_t i=0; i<b.size(); ++i) {
69    suite.out() << i << " " << b(i) << " " << stats[i].mean() << " "
70                << stats[i].standard_error() << "\n";
71    if (stats[i].standard_error() == 0.0) {
72      suite.xadd(false);
73      suite.err() << "error: standard error is 0.0\n";
74    }
75    else if (std::abs(stats[i].mean() - b(i)) > 5*stats[i].standard_error()) {
76      suite.xadd(false);
77      suite.err() << "error: average for param " << i << ": "
78                  << stats[i].mean() << "\n";
79    }
80  }
81
82  return suite.return_value();
83}
84
85
86void generate_data(const utility::Vector& b,
87                   utility::Matrix& X, utility::Vector& y)
88{
89  size_t n = 5000;
90  X.resize(n, b.size()-1);
91  y.resize(n);
92  random::Gaussian gauss;
93  random::Poisson poisson;
94  for (utility::Matrix::iterator it=X.begin(); it!=X.end(); ++it)
95    *it = gauss();
96
97  for (size_t i=0; i<n; ++i) {
98    double lnmu = b(0);
99    for (size_t j=1; j<b.size(); ++j)
100      lnmu += X(i, j-1) * b(j);
101
102    double mu = exp(lnmu);
103    y(i) = poisson(mu);
104  }
105}
106
107
108void analyse(const utility::Vector& b,
109             std::vector<statistics::Averager>& stats,
110             std::vector<statistics::Averager>& stats2)
111{
112  utility::Matrix X;
113  utility::Vector y;
114  generate_data(b, X, y);
115  theplu::yat::regression::Poisson poisson;
116  poisson.fit(X, y);
117  for (size_t i=0; i<poisson.fit_parameters().size(); ++i) {
118    stats[i].add(poisson.fit_parameters()(i));
119  }
120}
Note: See TracBrowser for help on using the repository browser.