test/poisson_mixture.cc

Code
Comments
Other
Rev Date Author Line
4180 09 Jun 22 peter 1 // $Id$
4180 09 Jun 22 peter 2
4180 09 Jun 22 peter 3 /*
4180 09 Jun 22 peter 4   Copyright (C) 2022 Peter Johansson
4180 09 Jun 22 peter 5
4180 09 Jun 22 peter 6   This file is part of the yat library, https://dev.thep.lu.se/yat
4180 09 Jun 22 peter 7
4180 09 Jun 22 peter 8   The yat library is free software; you can redistribute it and/or
4180 09 Jun 22 peter 9   modify it under the terms of the GNU General Public License as
4180 09 Jun 22 peter 10   published by the Free Software Foundation; either version 3 of the
4180 09 Jun 22 peter 11   License, or (at your option) any later version.
4180 09 Jun 22 peter 12
4180 09 Jun 22 peter 13   The yat library is distributed in the hope that it will be useful,
4180 09 Jun 22 peter 14   but WITHOUT ANY WARRANTY; without even the implied warranty of
4180 09 Jun 22 peter 15   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
4180 09 Jun 22 peter 16   General Public License for more details.
4180 09 Jun 22 peter 17
4180 09 Jun 22 peter 18   You should have received a copy of the GNU General Public License
4180 09 Jun 22 peter 19   along with yat. If not, see <https://www.gnu.org/licenses/>.
4180 09 Jun 22 peter 20 */
4180 09 Jun 22 peter 21
4180 09 Jun 22 peter 22 #include <config.h>
4180 09 Jun 22 peter 23
4180 09 Jun 22 peter 24 #include "Suite.h"
4180 09 Jun 22 peter 25
4180 09 Jun 22 peter 26 #include <yat/random/random.h>
4180 09 Jun 22 peter 27 #include <yat/statistics/PoissonMixture.h>
4180 09 Jun 22 peter 28
4184 30 Jun 22 peter 29 #include <cassert>
4184 30 Jun 22 peter 30
4180 09 Jun 22 peter 31 using namespace theplu::yat;
4180 09 Jun 22 peter 32
4180 09 Jun 22 peter 33 void test1(test::Suite& suite);
4180 09 Jun 22 peter 34
4180 09 Jun 22 peter 35 int main(int argc, char* argv[])
4180 09 Jun 22 peter 36 {
4180 09 Jun 22 peter 37   test::Suite suite(argc, argv);
4180 09 Jun 22 peter 38   try {
4180 09 Jun 22 peter 39     test1(suite);
4180 09 Jun 22 peter 40   }
4180 09 Jun 22 peter 41   catch (std::runtime_error& e) {
4180 09 Jun 22 peter 42     suite.err() << "error: " << e.what() << "\n";
4180 09 Jun 22 peter 43     suite.add(false);
4180 09 Jun 22 peter 44   }
4180 09 Jun 22 peter 45   return suite.return_value();
4180 09 Jun 22 peter 46 }
4180 09 Jun 22 peter 47
4180 09 Jun 22 peter 48
4180 09 Jun 22 peter 49 void test1(test::Suite& suite)
4180 09 Jun 22 peter 50 {
4180 09 Jun 22 peter 51   std::vector<double> means = {20, 50};
4180 09 Jun 22 peter 52   std::vector<double> tau(means.size(), 0.1);
4180 09 Jun 22 peter 53   statistics::Histogram hist(0, means.size(), means.size());
4180 09 Jun 22 peter 54   for (size_t i=0; i<tau.size(); ++i)
4180 09 Jun 22 peter 55     hist.add(i, tau[i]);
4180 09 Jun 22 peter 56   random::DiscreteGeneral rng(hist);
4180 09 Jun 22 peter 57   std::vector<random::Poisson> po;
4180 09 Jun 22 peter 58   po.reserve(means.size());
4180 09 Jun 22 peter 59   for (double m : means)
4180 09 Jun 22 peter 60     po.push_back(random::Poisson(m));
4180 09 Jun 22 peter 61
4180 09 Jun 22 peter 62   statistics::PoissonMixture model;
4180 09 Jun 22 peter 63   size_t N=10000;
4180 09 Jun 22 peter 64   std::vector<unsigned long int> data(N);
4180 09 Jun 22 peter 65   for (size_t i=0; i<10000; ++i) {
4180 09 Jun 22 peter 66     data[i] = po[rng()]();
4180 09 Jun 22 peter 67     model.add(data[i]);
4180 09 Jun 22 peter 68   }
4180 09 Jun 22 peter 69   model.fit(8);
4180 09 Jun 22 peter 70
4180 09 Jun 22 peter 71   suite.out() << "Mean\t\tTau\n";
4180 09 Jun 22 peter 72   suite.out() << "Inferred\tExpected\tInferred\tExpected\n";
4180 09 Jun 22 peter 73   for (size_t i=0; i<means.size(); ++i) {
4180 09 Jun 22 peter 74     suite.out() << model.mean(i) << "\t" << means[i] << "\t";
4180 09 Jun 22 peter 75     suite.out() << model.tau(i) << "\t" << tau[i] << "\n";
4180 09 Jun 22 peter 76   }
4180 09 Jun 22 peter 77   double log_L = model.logL();
4180 09 Jun 22 peter 78   suite.out() << "logL:  " << log_L << "\n";
4180 09 Jun 22 peter 79
4180 09 Jun 22 peter 80   double log_L0 = 0;
4180 09 Jun 22 peter 81   for (auto k : data) {
4180 09 Jun 22 peter 82     double L = 0;
4180 09 Jun 22 peter 83     for (size_t i=0; i<tau.size(); ++i)
4180 09 Jun 22 peter 84       L += tau[i] * gsl_ran_poisson_pdf(k, means[i]);
4180 09 Jun 22 peter 85     log_L0 += std::log(L);
4180 09 Jun 22 peter 86   }
4180 09 Jun 22 peter 87   suite.out() << "logL0: " << log_L0 << "\n";
4180 09 Jun 22 peter 88
4180 09 Jun 22 peter 89
4180 09 Jun 22 peter 90   if (log_L < log_L0) {
4180 09 Jun 22 peter 91     suite.add(false);
4180 09 Jun 22 peter 92     suite.err() << "logL for inferred model is smaller than for true model\n";
4180 09 Jun 22 peter 93   }
4180 09 Jun 22 peter 94
4180 09 Jun 22 peter 95 }

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