test/negative_binomial_extended_mixture.cc

Code
Comments
Other
Rev Date Author Line
4192 11 Aug 22 peter 1 // $Id$
4192 11 Aug 22 peter 2 //
4192 11 Aug 22 peter 3 // Copyright (C) 2022 Peter Johansson
4192 11 Aug 22 peter 4 //
4192 11 Aug 22 peter 5 // This program is free software; you can redistribute it and/or modify
4192 11 Aug 22 peter 6 // it under the terms of the GNU General Public License as published by
4192 11 Aug 22 peter 7 // the Free Software Foundation; either version 3 of the License, or
4192 11 Aug 22 peter 8 // (at your option) any later version.
4192 11 Aug 22 peter 9 //
4192 11 Aug 22 peter 10 // This program is distributed in the hope that it will be useful, but
4192 11 Aug 22 peter 11 // WITHOUT ANY WARRANTY; without even the implied warranty of
4192 11 Aug 22 peter 12 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
4192 11 Aug 22 peter 13 // General Public License for more details.
4192 11 Aug 22 peter 14 //
4192 11 Aug 22 peter 15 // You should have received a copy of the GNU General Public License
4192 11 Aug 22 peter 16 // along with this program. If not, see <https://www.gnu.org/licenses/>.
4192 11 Aug 22 peter 17
4192 11 Aug 22 peter 18 #include <config.h>
4192 11 Aug 22 peter 19
4192 11 Aug 22 peter 20 #include <cassert>
4192 11 Aug 22 peter 21
4192 11 Aug 22 peter 22 #include "Suite.h"
4192 11 Aug 22 peter 23
4192 11 Aug 22 peter 24 #include "yat/random/random.h"
4192 11 Aug 22 peter 25 #include "yat/statistics/NegativeBinomialExtendedMixture.h"
4192 11 Aug 22 peter 26
4192 11 Aug 22 peter 27 using namespace theplu::yat;
4192 11 Aug 22 peter 28
4192 11 Aug 22 peter 29 void test1(test::Suite& suite);
4192 11 Aug 22 peter 30
4192 11 Aug 22 peter 31 int main(int argc, char* argv[])
4192 11 Aug 22 peter 32 {
4192 11 Aug 22 peter 33   test::Suite suite(argc, argv);
4192 11 Aug 22 peter 34   try {
4192 11 Aug 22 peter 35     test1(suite);
4192 11 Aug 22 peter 36   }
4192 11 Aug 22 peter 37   catch (std::runtime_error& e) {
4192 11 Aug 22 peter 38     suite.err() << "error: " << e.what() << "\n";
4192 11 Aug 22 peter 39     return EXIT_FAILURE;
4192 11 Aug 22 peter 40   }
4192 11 Aug 22 peter 41   return suite.return_value();
4192 11 Aug 22 peter 42 }
4192 11 Aug 22 peter 43
4192 11 Aug 22 peter 44
4192 11 Aug 22 peter 45 void test1(test::Suite& suite)
4192 11 Aug 22 peter 46 {
4192 11 Aug 22 peter 47   double alpha = 0.3;
4192 11 Aug 22 peter 48   utility::Vector m(2);
4192 11 Aug 22 peter 49   m(0) = 100;
4192 11 Aug 22 peter 50   m(1) = 400;
4192 11 Aug 22 peter 51   utility::Vector a(2);
4192 11 Aug 22 peter 52   a(0) = 1.01;
4192 11 Aug 22 peter 53   a(1) = 2;
4192 11 Aug 22 peter 54
4192 11 Aug 22 peter 55   utility::Vector r(2);
4192 11 Aug 22 peter 56   utility::Vector p(2);
4192 11 Aug 22 peter 57   for (size_t i=0; i<2; ++i) {
4192 11 Aug 22 peter 58     r(i) = m(i)/(a(i)-1);
4192 11 Aug 22 peter 59     p(i) = 1 - 1/a(i);
4192 11 Aug 22 peter 60   }
4192 11 Aug 22 peter 61
4192 11 Aug 22 peter 62   for (size_t i=0; i<m.size(); ++i) {
4192 11 Aug 22 peter 63     suite.out() << i << "\t";
4192 11 Aug 22 peter 64     suite.out() << "m: " << m(i) << " ";
4192 11 Aug 22 peter 65     suite.out() << "alpha: " << a(i) << "\n";
4192 11 Aug 22 peter 66   }
4192 11 Aug 22 peter 67
4192 11 Aug 22 peter 68   statistics::NegativeBinomialExtendedMixture nbem;
4192 11 Aug 22 peter 69   random::Exponential exponential;
4192 11 Aug 22 peter 70   random::ContinuousUniform urng;
4192 11 Aug 22 peter 71   std::vector<std::pair<unsigned long int, double>> data;
4192 11 Aug 22 peter 72   for (size_t i=0; i<10000; ++i) {
4192 11 Aug 22 peter 73     double y = exponential();
4192 11 Aug 22 peter 74     //y = 1.0;
4192 11 Aug 22 peter 75     size_t idx = 0;
4192 11 Aug 22 peter 76     if (urng() < alpha)
4192 11 Aug 22 peter 77       idx = 0;
4192 11 Aug 22 peter 78     else
4192 11 Aug 22 peter 79       idx = 1;
4192 11 Aug 22 peter 80
4192 11 Aug 22 peter 81     random::NegativeBinomial rng(1.0 - p(idx), y*r(idx));
4192 11 Aug 22 peter 82     double x = rng();
4192 11 Aug 22 peter 83     data.push_back(std::make_pair(x, y));
4192 11 Aug 22 peter 84     nbem.add(x, y);
4192 11 Aug 22 peter 85   }
4192 11 Aug 22 peter 86   nbem.fit(2);
4192 11 Aug 22 peter 87
4192 11 Aug 22 peter 88   suite.out() << "tau:  " << nbem.tau() << "\n";
4192 11 Aug 22 peter 89   if (std::abs(nbem.tau()(0) - alpha) > 0.01) {
4192 11 Aug 22 peter 90     suite.add(false);
4192 11 Aug 22 peter 91     suite.err() << "error: incorrect: tau: " << nbem.tau() << "\n";
4192 11 Aug 22 peter 92   }
4192 11 Aug 22 peter 93   suite.out() << "m:    " << nbem.m() << "\n";
4192 11 Aug 22 peter 94   for (size_t i=0; i<2; ++i) {
4192 11 Aug 22 peter 95     if (std::abs(nbem.m()(i) - m(i)) > 1.0) {
4192 11 Aug 22 peter 96       suite.err() << "incorrect m: " << nbem.m()(i) << "\n";
4192 11 Aug 22 peter 97       suite.add(false);
4192 11 Aug 22 peter 98     }
4192 11 Aug 22 peter 99   }
4192 11 Aug 22 peter 100   suite.out() << "alpha:" << nbem.alpha() << "\n";
4192 11 Aug 22 peter 101
4192 11 Aug 22 peter 102   double logL = nbem.logL();
4192 11 Aug 22 peter 103   suite.out() << "logL:  " << logL << "\n";
4192 11 Aug 22 peter 104
4192 11 Aug 22 peter 105   double logL0 = 0;
4192 11 Aug 22 peter 106   for (auto x : data) {
4192 11 Aug 22 peter 107     unsigned long int k = x.first;
4192 11 Aug 22 peter 108     double y = x.second;
4192 11 Aug 22 peter 109
4192 11 Aug 22 peter 110     double pdf = 0;
4192 11 Aug 22 peter 111     for (size_t i=0; i<2; ++i) {
4192 11 Aug 22 peter 112       double tau = i ? 1-alpha : alpha;
4192 11 Aug 22 peter 113       pdf += tau * gsl_ran_negative_binomial_pdf(k, 1.0-p(i), y*r(i));
4192 11 Aug 22 peter 114     }
4192 11 Aug 22 peter 115     logL0 += std::log(pdf);
4192 11 Aug 22 peter 116   }
4192 11 Aug 22 peter 117   suite.out() << "logL0: " << logL0 << "\n";
4192 11 Aug 22 peter 118
4192 11 Aug 22 peter 119   if (logL < logL0) {
4192 11 Aug 22 peter 120     suite.add(false);
4192 11 Aug 22 peter 121     suite.err() << "logL for inferred model is smaller than for true model\n";
4192 11 Aug 22 peter 122   }
4192 11 Aug 22 peter 123
4192 11 Aug 22 peter 124 }

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