test/negative_binomial.cc

Code
Comments
Other
Rev Date Author Line
3659 13 Jul 17 peter 1 // $Id$
3659 13 Jul 17 peter 2
3659 13 Jul 17 peter 3 /*
3999 08 Oct 20 peter 4   Copyright (C) 2017, 2019, 2020 Peter Johansson
3659 13 Jul 17 peter 5
3659 13 Jul 17 peter 6   This file is part of the yat library, http://dev.thep.lu.se/yat
3659 13 Jul 17 peter 7
3659 13 Jul 17 peter 8   The yat library is free software; you can redistribute it and/or
3659 13 Jul 17 peter 9   modify it under the terms of the GNU General Public License as
3659 13 Jul 17 peter 10   published by the Free Software Foundation; either version 3 of the
3659 13 Jul 17 peter 11   License, or (at your option) any later version.
3659 13 Jul 17 peter 12
3659 13 Jul 17 peter 13   The yat library is distributed in the hope that it will be useful,
3659 13 Jul 17 peter 14   but WITHOUT ANY WARRANTY; without even the implied warranty of
3659 13 Jul 17 peter 15   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
3659 13 Jul 17 peter 16   General Public License for more details.
3659 13 Jul 17 peter 17
3659 13 Jul 17 peter 18   You should have received a copy of the GNU General Public License
3659 13 Jul 17 peter 19   along with yat. If not, see <http://www.gnu.org/licenses/>.
3659 13 Jul 17 peter 20 */
3659 13 Jul 17 peter 21
3659 13 Jul 17 peter 22 #include <config.h>
3659 13 Jul 17 peter 23
3659 13 Jul 17 peter 24 #include "Suite.h"
3659 13 Jul 17 peter 25
3659 13 Jul 17 peter 26 #include "yat/regression/NegativeBinomial.h"
3659 13 Jul 17 peter 27
3659 13 Jul 17 peter 28 #include "yat/random/random.h"
3659 13 Jul 17 peter 29 #include "yat/statistics/Averager.h"
3659 13 Jul 17 peter 30
3659 13 Jul 17 peter 31 #include "yat/utility/Matrix.h"
3659 13 Jul 17 peter 32 #include "yat/utility/Vector.h"
3659 13 Jul 17 peter 33
3788 29 Jan 19 peter 34 #include <cmath>
3659 13 Jul 17 peter 35 #include <vector>
3659 13 Jul 17 peter 36
3659 13 Jul 17 peter 37 using namespace theplu::yat;
3659 13 Jul 17 peter 38
3659 13 Jul 17 peter 39 void analyse(const utility::Vector& b,
3659 13 Jul 17 peter 40              std::vector<statistics::Averager>& stats,
3659 13 Jul 17 peter 41              std::vector<statistics::Averager>& stats2,
3663 20 Jul 17 peter 42              statistics::Averager& stats3,
3663 20 Jul 17 peter 43              double alpha, const utility::Matrix&);
3659 13 Jul 17 peter 44
3663 20 Jul 17 peter 45 void generate_X(const utility::Vector& b, utility::Matrix& X, double alpha);
3663 20 Jul 17 peter 46
3980 24 Aug 20 peter 47
3980 24 Aug 20 peter 48 void some_function(regression::Multivariate&)
3980 24 Aug 20 peter 49 {
3980 24 Aug 20 peter 50 }
3980 24 Aug 20 peter 51
3980 24 Aug 20 peter 52
3659 13 Jul 17 peter 53 int main(int argc, char* argv[])
3659 13 Jul 17 peter 54 {
3659 13 Jul 17 peter 55   test::Suite suite(argc, argv);
3659 13 Jul 17 peter 56
3659 13 Jul 17 peter 57   double alpha = 2.0;
3659 13 Jul 17 peter 58   utility::Vector b(4);
3659 13 Jul 17 peter 59   b(0) = 0.5;
3659 13 Jul 17 peter 60   b(1) = 2.0;
3659 13 Jul 17 peter 61   b(2) = 0.75;
3659 13 Jul 17 peter 62   b(3) = -1.25;
3659 13 Jul 17 peter 63   std::vector<statistics::Averager> stats(b.size());
3659 13 Jul 17 peter 64   std::vector<statistics::Averager> stats2(b.size());
3663 20 Jul 17 peter 65   statistics::Averager stats3;
3663 20 Jul 17 peter 66   utility::Matrix X;
3663 20 Jul 17 peter 67   generate_X(b, X, alpha);
3659 13 Jul 17 peter 68   for (size_t i=0; i<100; ++i)
3663 20 Jul 17 peter 69     analyse(b, stats, stats2, stats3, alpha, X);
3659 13 Jul 17 peter 70
3663 20 Jul 17 peter 71   suite.out() << "\ttrue value\tmean\tz\n";
3663 20 Jul 17 peter 72   for (size_t i=0; i<b.size(); ++i) {
3663 20 Jul 17 peter 73     double z = (stats[i].mean()-b(i)) / stats[i].standard_error();
3663 20 Jul 17 peter 74     suite.out() << i << "\t"
3663 20 Jul 17 peter 75                 << b(i) << "\t"
3663 20 Jul 17 peter 76                 << stats[i].mean() << "\t"
3663 20 Jul 17 peter 77                 << z << "\n";
3663 20 Jul 17 peter 78     if (stats[i].standard_error() == 0.0) {
3663 20 Jul 17 peter 79       suite.add(false);
3663 20 Jul 17 peter 80       suite.err() << "error: standard error is 0.0\n";
3663 20 Jul 17 peter 81     }
3663 20 Jul 17 peter 82     else if (std::abs(z) > 5) {
3663 20 Jul 17 peter 83       suite.add(false);
3663 20 Jul 17 peter 84       suite.err() << "error: average for param " << i << ": "
3663 20 Jul 17 peter 85                   << stats[i].mean() << "\n";
3663 20 Jul 17 peter 86     }
3663 20 Jul 17 peter 87   }
3663 20 Jul 17 peter 88
3663 20 Jul 17 peter 89   suite.out() << "\ntest covariance\n";
3663 20 Jul 17 peter 90   suite.out() << "i\t<inferred>\tobserved\tdelta\trelative error\n";
3663 20 Jul 17 peter 91   for (size_t i=0; i<b.size(); ++i) {
3663 20 Jul 17 peter 92     // comparing uncertainty observed (stats2) and uncertainty
3663 20 Jul 17 peter 93     // estimated (stats)
3663 20 Jul 17 peter 94
3663 20 Jul 17 peter 95     double delta = stats2[i].mean() - stats[i].std();
3663 20 Jul 17 peter 96     double relative_error = delta / stats[i].std();
3663 20 Jul 17 peter 97
3663 20 Jul 17 peter 98     suite.out() << i << "\t" << stats2[i].mean() << "\t"
3663 20 Jul 17 peter 99                 << stats[i].std() << "\t"
3663 20 Jul 17 peter 100                 << delta << "\t"
3663 20 Jul 17 peter 101                 << relative_error << "\n";
3663 20 Jul 17 peter 102
3663 20 Jul 17 peter 103     // covariance matrix is only an approximation, so only check that
3663 20 Jul 17 peter 104     // estimation is roughly correct (within 40%)
3663 20 Jul 17 peter 105     if (std::abs(relative_error) > 0.4) {
3663 20 Jul 17 peter 106       suite.add(false);
3663 20 Jul 17 peter 107       suite.err() << "error: " << relative_error << " too large\n";
3663 20 Jul 17 peter 108     }
3663 20 Jul 17 peter 109   }
3663 20 Jul 17 peter 110
3663 20 Jul 17 peter 111   suite.out() << "\nalpha: " << stats3.mean() << "\t"
3663 20 Jul 17 peter 112               << stats3.std() << "\t"
3663 20 Jul 17 peter 113               << alpha
3663 20 Jul 17 peter 114               << "\n";
3663 20 Jul 17 peter 115
3788 29 Jan 19 peter 116   if (std::fabs(stats3.mean() - alpha) > stats3.std()) {
3663 20 Jul 17 peter 117     suite.add(false);
3663 20 Jul 17 peter 118     suite.err() << "incorrect alpha\n";
3663 20 Jul 17 peter 119   }
3663 20 Jul 17 peter 120
3659 13 Jul 17 peter 121   return suite.return_value();
3659 13 Jul 17 peter 122 }
3659 13 Jul 17 peter 123
3663 20 Jul 17 peter 124
3663 20 Jul 17 peter 125 void generate_X(const utility::Vector& b, utility::Matrix& X, double alpha)
3659 13 Jul 17 peter 126 {
3659 13 Jul 17 peter 127   size_t n = 5000;
3659 13 Jul 17 peter 128   X.resize(n, b.size()-1);
3659 13 Jul 17 peter 129   random::Gaussian gauss;
3659 13 Jul 17 peter 130   for (utility::Matrix::iterator it=X.begin(); it!=X.end(); ++it)
3659 13 Jul 17 peter 131     *it = gauss();
3663 20 Jul 17 peter 132 }
3659 13 Jul 17 peter 133
3663 20 Jul 17 peter 134
3663 20 Jul 17 peter 135 void generate_y(const utility::Vector& b, const utility::Matrix& X,
3663 20 Jul 17 peter 136                 utility::Vector& y, double alpha)
3663 20 Jul 17 peter 137 {
3663 20 Jul 17 peter 138   size_t n = X.rows();
3663 20 Jul 17 peter 139   y.resize(n);
3663 20 Jul 17 peter 140
3659 13 Jul 17 peter 141   for (size_t i=0; i<n; ++i) {
3659 13 Jul 17 peter 142     double lnmu = b(0);
3659 13 Jul 17 peter 143     for (size_t j=1; j<b.size(); ++j)
3659 13 Jul 17 peter 144       lnmu += X(i, j-1) * b(j);
3659 13 Jul 17 peter 145
3659 13 Jul 17 peter 146     double mu = exp(lnmu);
3659 13 Jul 17 peter 147
3659 13 Jul 17 peter 148     // alpha = 1.0 / (1-p) -> p = 1 - 1/alpha
3659 13 Jul 17 peter 149     double p = 1.0 - 1.0 / alpha;
3663 20 Jul 17 peter 150     // mu = pr / (1-p) -> r = mu (1-p) / p
3659 13 Jul 17 peter 151     double r = mu * (1-p) / p;
3659 13 Jul 17 peter 152
3659 13 Jul 17 peter 153
3659 13 Jul 17 peter 154     // GSL turns things upside down and defines as follows
3659 13 Jul 17 peter 155     //
3659 13 Jul 17 peter 156     // This function returns "the number of failures occurring before
3659 13 Jul 17 peter 157     // n successes in independent trials with probability p of
3659 13 Jul 17 peter 158     // success"
3659 13 Jul 17 peter 159     //
3659 13 Jul 17 peter 160     //   gsl_ran_negative_binomial(RNG->rng(), p, n);
3659 13 Jul 17 peter 161     //
3659 13 Jul 17 peter 162     // thus
3659 13 Jul 17 peter 163     //
3659 13 Jul 17 peter 164     //   gsl_ran_negative_binomial(RNG->rng(), 1-p, r);
3659 13 Jul 17 peter 165     //
3659 13 Jul 17 peter 166     // returns number of failures occurring before r sucsesses in
3659 13 Jul 17 peter 167     // independent trials with success probability 1-p, or if we swap
3659 13 Jul 17 peter 168     // success and failure: number of sucesses before r failures in
3659 13 Jul 17 peter 169     // independent trials with success probability p.
3659 13 Jul 17 peter 170
3659 13 Jul 17 peter 171     y(i) = gsl_ran_negative_binomial(random::RNG::instance()->rng(), 1-p, r);
3659 13 Jul 17 peter 172   }
3659 13 Jul 17 peter 173 }
3659 13 Jul 17 peter 174
3659 13 Jul 17 peter 175
3659 13 Jul 17 peter 176 void analyse(const utility::Vector& b,
3659 13 Jul 17 peter 177              std::vector<statistics::Averager>& stats,
3659 13 Jul 17 peter 178              std::vector<statistics::Averager>& stats2,
3663 20 Jul 17 peter 179              statistics::Averager& stats3,
3663 20 Jul 17 peter 180              double alpha, const utility::Matrix& X)
3659 13 Jul 17 peter 181 {
3659 13 Jul 17 peter 182   utility::Vector y;
3663 20 Jul 17 peter 183   generate_y(b, X, y, alpha);
3663 20 Jul 17 peter 184   regression::NegativeBinomial regress;
3663 20 Jul 17 peter 185   regress.fit(X, y);
3663 20 Jul 17 peter 186   for (size_t i=0; i<regress.fit_parameters().size(); ++i) {
3663 20 Jul 17 peter 187     stats[i].add(regress.fit_parameters()(i));
3663 20 Jul 17 peter 188     stats2[i].add(std::sqrt(regress.covariance()(i, i)));
3663 20 Jul 17 peter 189   }
3663 20 Jul 17 peter 190   stats3.add(regress.alpha());
3980 24 Aug 20 peter 191   some_function(regress);
3659 13 Jul 17 peter 192 }

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