test/poisson.cc

Code
Comments
Other
Rev Date Author Line
3615 06 Feb 17 peter 1 // $Id$
3615 06 Feb 17 peter 2
3615 06 Feb 17 peter 3 /*
3648 02 Jun 17 peter 4   Copyright (C) 2017 Peter Johansson
3615 06 Feb 17 peter 5
3615 06 Feb 17 peter 6   This file is part of the yat library, http://dev.thep.lu.se/yat
3615 06 Feb 17 peter 7
3615 06 Feb 17 peter 8   The yat library is free software; you can redistribute it and/or
3615 06 Feb 17 peter 9   modify it under the terms of the GNU General Public License as
3615 06 Feb 17 peter 10   published by the Free Software Foundation; either version 3 of the
3615 06 Feb 17 peter 11   License, or (at your option) any later version.
3615 06 Feb 17 peter 12
3615 06 Feb 17 peter 13   The yat library is distributed in the hope that it will be useful,
3615 06 Feb 17 peter 14   but WITHOUT ANY WARRANTY; without even the implied warranty of
3615 06 Feb 17 peter 15   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
3615 06 Feb 17 peter 16   General Public License for more details.
3615 06 Feb 17 peter 17
3615 06 Feb 17 peter 18   You should have received a copy of the GNU General Public License
3615 06 Feb 17 peter 19   along with yat. If not, see <http://www.gnu.org/licenses/>.
3615 06 Feb 17 peter 20 */
3615 06 Feb 17 peter 21
3615 06 Feb 17 peter 22 #include <config.h>
3615 06 Feb 17 peter 23
3615 06 Feb 17 peter 24 #include "Suite.h"
3615 06 Feb 17 peter 25
3615 06 Feb 17 peter 26 #include "yat/regression/Poisson.h"
3648 02 Jun 17 peter 27 #include "yat/random/random.h"
3648 02 Jun 17 peter 28 #include "yat/statistics/Averager.h"
3615 06 Feb 17 peter 29 #include "yat/utility/Matrix.h"
3615 06 Feb 17 peter 30 #include "yat/utility/Vector.h"
3615 06 Feb 17 peter 31
3648 02 Jun 17 peter 32 #include <vector>
3648 02 Jun 17 peter 33
3615 06 Feb 17 peter 34 using namespace theplu::yat;
3615 06 Feb 17 peter 35
3662 19 Jul 17 peter 36 void analyse(const utility::Vector& b, const utility::Matrix& X,
3648 02 Jun 17 peter 37              std::vector<statistics::Averager>& stats,
3648 02 Jun 17 peter 38              std::vector<statistics::Averager>& stats2);
3648 02 Jun 17 peter 39
3662 19 Jul 17 peter 40 void generate_X(const utility::Vector& b, utility::Matrix& X);
3662 19 Jul 17 peter 41
3615 06 Feb 17 peter 42 int main(int argc, char* argv[])
3615 06 Feb 17 peter 43 {
3615 06 Feb 17 peter 44   test::Suite suite(argc, argv);
3615 06 Feb 17 peter 45
3615 06 Feb 17 peter 46   regression::Poisson model;
3662 19 Jul 17 peter 47   size_t n = 1000;
3615 06 Feb 17 peter 48   size_t p = 4;
3615 06 Feb 17 peter 49   utility::Vector y(n);
3615 06 Feb 17 peter 50   utility::Matrix X(n, p);
3615 06 Feb 17 peter 51   model.fit(X, y);
3615 06 Feb 17 peter 52   if (model.fit_parameters().size() != p+1) {
3615 06 Feb 17 peter 53     suite.add(false);
3615 06 Feb 17 peter 54     suite.err() << "error: size of fit_parameters: "
3615 06 Feb 17 peter 55                 << model.fit_parameters().size()
3615 06 Feb 17 peter 56                 << " expected " << p+1 << "\n";
3615 06 Feb 17 peter 57   }
3615 06 Feb 17 peter 58   model.predict(X.row_const_view(0));
3615 06 Feb 17 peter 59
3648 02 Jun 17 peter 60   utility::Vector b(4);
3648 02 Jun 17 peter 61   b(0) = 0.5;
3648 02 Jun 17 peter 62   b(1) = 2.0;
3648 02 Jun 17 peter 63   b(2) = 0.75;
3648 02 Jun 17 peter 64   b(3) = -1.25;
3662 19 Jul 17 peter 65   generate_X(b, X);
3648 02 Jun 17 peter 66   std::vector<statistics::Averager> stats(b.size());
3648 02 Jun 17 peter 67   std::vector<statistics::Averager> stats2(b.size());
3648 02 Jun 17 peter 68   for (size_t i=0; i<100; ++i)
3662 19 Jul 17 peter 69     analyse(b, X, stats, stats2);
3648 02 Jun 17 peter 70
3648 02 Jun 17 peter 71   for (size_t i=0; i<b.size(); ++i) {
3648 02 Jun 17 peter 72     suite.out() << i << " " << b(i) << " " << stats[i].mean() << " "
3648 02 Jun 17 peter 73                 << stats[i].standard_error() << "\n";
3648 02 Jun 17 peter 74     if (stats[i].standard_error() == 0.0) {
3658 13 Jul 17 peter 75       suite.add(false);
3648 02 Jun 17 peter 76       suite.err() << "error: standard error is 0.0\n";
3648 02 Jun 17 peter 77     }
3648 02 Jun 17 peter 78     else if (std::abs(stats[i].mean() - b(i)) > 5*stats[i].standard_error()) {
3658 13 Jul 17 peter 79       suite.add(false);
3648 02 Jun 17 peter 80       suite.err() << "error: average for param " << i << ": "
3648 02 Jun 17 peter 81                   << stats[i].mean() << "\n";
3648 02 Jun 17 peter 82     }
3648 02 Jun 17 peter 83   }
3648 02 Jun 17 peter 84
3662 19 Jul 17 peter 85   suite.out() << "\ntest covariance\n";
3662 19 Jul 17 peter 86   suite.out() << "i\t<inferred>\tobserved\tdelta\trelative error\n";
3662 19 Jul 17 peter 87   for (size_t i=0; i<b.size(); ++i) {
3662 19 Jul 17 peter 88     // comparing uncertainty observed (stats2) and uncertainty
3662 19 Jul 17 peter 89     // estimated (stats)
3662 19 Jul 17 peter 90
3662 19 Jul 17 peter 91     double delta = stats2[i].mean() - stats[i].std();
3662 19 Jul 17 peter 92     double relative_error = delta / stats[i].std();
3662 19 Jul 17 peter 93
3662 19 Jul 17 peter 94     suite.out() << i << "\t" << stats2[i].mean() << "\t"
3662 19 Jul 17 peter 95                 << stats[i].std() << "\t"
3662 19 Jul 17 peter 96                 << delta << "\t"
3662 19 Jul 17 peter 97                 << relative_error << "\n";
3662 19 Jul 17 peter 98
3662 19 Jul 17 peter 99     // covariance matrix is only an approximation, so only check that
3662 19 Jul 17 peter 100     // estimation is roughly correct (within 20%)
3662 19 Jul 17 peter 101     if (std::abs(relative_error) > 0.2) {
3662 19 Jul 17 peter 102       suite.add(false);
3662 19 Jul 17 peter 103       suite.err() << "error: " << relative_error << " too large\n";
3662 19 Jul 17 peter 104     }
3662 19 Jul 17 peter 105   }
3662 19 Jul 17 peter 106
3615 06 Feb 17 peter 107   return suite.return_value();
3615 06 Feb 17 peter 108 }
3648 02 Jun 17 peter 109
3648 02 Jun 17 peter 110
3662 19 Jul 17 peter 111 void generate_X(const utility::Vector& b, utility::Matrix& X)
3648 02 Jun 17 peter 112 {
3648 02 Jun 17 peter 113   size_t n = 5000;
3648 02 Jun 17 peter 114   X.resize(n, b.size()-1);
3648 02 Jun 17 peter 115   random::Gaussian gauss;
3648 02 Jun 17 peter 116   for (utility::Matrix::iterator it=X.begin(); it!=X.end(); ++it)
3648 02 Jun 17 peter 117     *it = gauss();
3662 19 Jul 17 peter 118 }
3648 02 Jun 17 peter 119
3662 19 Jul 17 peter 120
3662 19 Jul 17 peter 121 void generate_y(const utility::Vector& b,
3662 19 Jul 17 peter 122                 const utility::Matrix& X, utility::Vector& y)
3662 19 Jul 17 peter 123 {
3662 19 Jul 17 peter 124   size_t n = X.rows();
3662 19 Jul 17 peter 125   y.resize(n);
3662 19 Jul 17 peter 126   random::Poisson poisson;
3662 19 Jul 17 peter 127
3648 02 Jun 17 peter 128   for (size_t i=0; i<n; ++i) {
3648 02 Jun 17 peter 129     double lnmu = b(0);
3648 02 Jun 17 peter 130     for (size_t j=1; j<b.size(); ++j)
3648 02 Jun 17 peter 131       lnmu += X(i, j-1) * b(j);
3648 02 Jun 17 peter 132
3648 02 Jun 17 peter 133     double mu = exp(lnmu);
3648 02 Jun 17 peter 134     y(i) = poisson(mu);
3648 02 Jun 17 peter 135   }
3648 02 Jun 17 peter 136 }
3648 02 Jun 17 peter 137
3648 02 Jun 17 peter 138
3662 19 Jul 17 peter 139 void analyse(const utility::Vector& b, const utility::Matrix& X,
3648 02 Jun 17 peter 140              std::vector<statistics::Averager>& stats,
3648 02 Jun 17 peter 141              std::vector<statistics::Averager>& stats2)
3648 02 Jun 17 peter 142 {
3648 02 Jun 17 peter 143   utility::Vector y;
3662 19 Jul 17 peter 144   generate_y(b, X, y);
3648 02 Jun 17 peter 145   theplu::yat::regression::Poisson poisson;
3648 02 Jun 17 peter 146   poisson.fit(X, y);
3648 02 Jun 17 peter 147   for (size_t i=0; i<poisson.fit_parameters().size(); ++i) {
3648 02 Jun 17 peter 148     stats[i].add(poisson.fit_parameters()(i));
3662 19 Jul 17 peter 149     stats2[i].add(std::sqrt(poisson.covariance()(i, i)));
3648 02 Jun 17 peter 150   }
3648 02 Jun 17 peter 151 }

Generated on Sat Oct 14 02:21:22 2023 (UTC) by svndigest 0.11.2pre (dev build 1672)