test/negative_binomial_mixture.cc

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

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