test/likelihood.cc

Code
Comments
Other
Rev Date Author Line
4289 02 Feb 23 peter 1 // $Id$
4289 02 Feb 23 peter 2
4289 02 Feb 23 peter 3 /*
4289 02 Feb 23 peter 4   Copyright (C) 2023 Peter Johansson
4289 02 Feb 23 peter 5
4289 02 Feb 23 peter 6   This file is part of the yat library, https://dev.thep.lu.se/trac/yat
4289 02 Feb 23 peter 7
4289 02 Feb 23 peter 8   The yat library is free software; you can redistribute it and/or
4289 02 Feb 23 peter 9   modify it under the terms of the GNU General Public License as
4289 02 Feb 23 peter 10   published by the Free Software Foundation; either version 3 of the
4289 02 Feb 23 peter 11   License, or (at your option) any later version.
4289 02 Feb 23 peter 12
4289 02 Feb 23 peter 13   The yat library is distributed in the hope that it will be useful,
4289 02 Feb 23 peter 14   but WITHOUT ANY WARRANTY; without even the implied warranty of
4289 02 Feb 23 peter 15   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
4289 02 Feb 23 peter 16   General Public License for more details.
4289 02 Feb 23 peter 17
4289 02 Feb 23 peter 18   You should have received a copy of the GNU General Public License
4289 02 Feb 23 peter 19   along with yat. If not, see <https://www.gnu.org/licenses/>.
4289 02 Feb 23 peter 20 */
4289 02 Feb 23 peter 21
4289 02 Feb 23 peter 22 #include <config.h>
4289 02 Feb 23 peter 23
4289 02 Feb 23 peter 24 #include "Suite.h"
4289 02 Feb 23 peter 25
4289 02 Feb 23 peter 26 #include "yat/statistics/Likelihood.h"
4289 02 Feb 23 peter 27
4289 02 Feb 23 peter 28 #include <limits>
4289 02 Feb 23 peter 29 #include <vector>
4289 02 Feb 23 peter 30
4289 02 Feb 23 peter 31 using namespace theplu::yat;
4289 02 Feb 23 peter 32
4289 02 Feb 23 peter 33 int main(int argc, char* argv[])
4289 02 Feb 23 peter 34 {
4289 02 Feb 23 peter 35   test::Suite suite(argc, argv);
4289 02 Feb 23 peter 36   // ln(sum(exp(x))
4289 02 Feb 23 peter 37
4289 02 Feb 23 peter 38   statistics::Likelihood L;
4289 02 Feb 23 peter 39   for (double logL : { 0.01, 2.42, 1024.13, -12.34, 1020.13 })
4289 02 Feb 23 peter 40     L += statistics::exp(logL);
4289 02 Feb 23 peter 41
4289 02 Feb 23 peter 42   suite.out() << "log L: " << L.log() << "\n";
4289 02 Feb 23 peter 43   // we expect the sum to be slightly larger than max logL
4289 02 Feb 23 peter 44   if (L.log() <= 1024.13) {
4289 02 Feb 23 peter 45     suite.add(false);
4289 02 Feb 23 peter 46     suite.err() << "error: log L too small\n";
4289 02 Feb 23 peter 47   }
4289 02 Feb 23 peter 48   if (L.log() > 1024.16) {
4289 02 Feb 23 peter 49     suite.add(false);
4289 02 Feb 23 peter 50     suite.err() << "error: log L too large\n";
4289 02 Feb 23 peter 51   }
4289 02 Feb 23 peter 52
4289 02 Feb 23 peter 53   using statistics::Likelihood;
4289 02 Feb 23 peter 54   if (!suite.equal(log(Likelihood(1.23)),
4289 02 Feb 23 peter 55                    log(Likelihood(1.0)+Likelihood(0.23)))) {
4289 02 Feb 23 peter 56     suite.add(false);
4289 02 Feb 23 peter 57     suite.err() << "operator+ failed\n";
4289 02 Feb 23 peter 58   }
4289 02 Feb 23 peter 59   if (!suite.equal(log(2*Likelihood(1.23)), log(Likelihood(2*1.23)))) {
4289 02 Feb 23 peter 60     suite.add(false);
4289 02 Feb 23 peter 61     suite.err() << "operator* failed\n";
4289 02 Feb 23 peter 62   }
4289 02 Feb 23 peter 63
4289 02 Feb 23 peter 64   statistics::Likelihood L1(1.0);
4289 02 Feb 23 peter 65   L1 *= 10.2;
4289 02 Feb 23 peter 66   L1 /= 2.0;
4289 02 Feb 23 peter 67   if (!suite.equal(L1.get(), 5.1)) {
4289 02 Feb 23 peter 68     suite.add(false);
4289 02 Feb 23 peter 69     suite.err() << "operator*= or /= failed\n";
4289 02 Feb 23 peter 70   }
4289 02 Feb 23 peter 71
4289 02 Feb 23 peter 72   statistics::Likelihood L2;
4289 02 Feb 23 peter 73   if (!suite.add(L2.get() == 0.0)) {
4289 02 Feb 23 peter 74     suite.err() << "get() not zero\n";
4289 02 Feb 23 peter 75     suite.err() << "get(): " << L2.get() << " log(): " << L2.log() << "\n";
4289 02 Feb 23 peter 76   }
4289 02 Feb 23 peter 77   statistics::Likelihood L3(0.0);
4289 02 Feb 23 peter 78   if (!suite.add(L2 == L3)) {
4289 02 Feb 23 peter 79     suite.err() << "Likelihood() does not equal Likelihood(0.0)\n";
4289 02 Feb 23 peter 80     suite.err() << "logL2: " << L2.log() << " logL3: " << L3.log() << "\n";
4289 02 Feb 23 peter 81   }
4289 02 Feb 23 peter 82
4289 02 Feb 23 peter 83   std::vector<statistics::Likelihood> likelihoods;
4289 02 Feb 23 peter 84   likelihoods.push_back(Likelihood(12));
4289 02 Feb 23 peter 85   likelihoods.push_back(Likelihood(0));
4289 02 Feb 23 peter 86   likelihoods.push_back(Likelihood(1));
4289 02 Feb 23 peter 87   likelihoods.push_back(Likelihood(1000));
4289 02 Feb 23 peter 88   std::sort(likelihoods.begin(), likelihoods.end());
4289 02 Feb 23 peter 89   suite.add(suite.equal(likelihoods[0].get(), 0.0));
4289 02 Feb 23 peter 90   suite.add(suite.equal(likelihoods[1].get(), 1.0));
4289 02 Feb 23 peter 91   suite.add(suite.equal(likelihoods[2].get(), 12));
4378 11 Oct 23 peter 92   suite.add(suite.equal(likelihoods[3].get(), 1000, 22));
4289 02 Feb 23 peter 93
4289 02 Feb 23 peter 94   return suite.return_value();
4289 02 Feb 23 peter 95 }