test/chi2.cc

Code
Comments
Other
Rev Date Author Line
3723 05 Feb 18 peter 1 // $Id$
3723 05 Feb 18 peter 2
3723 05 Feb 18 peter 3 /*
3723 05 Feb 18 peter 4   Copyright (C) 2018 Peter Johansson
3723 05 Feb 18 peter 5
3723 05 Feb 18 peter 6   This file is part of the yat library, http://dev.thep.lu.se/yat
3723 05 Feb 18 peter 7
3723 05 Feb 18 peter 8   The yat library is free software; you can redistribute it and/or
3723 05 Feb 18 peter 9   modify it under the terms of the GNU General Public License as
3723 05 Feb 18 peter 10   published by the Free Software Foundation; either version 3 of the
3723 05 Feb 18 peter 11   License, or (at your option) any later version.
3723 05 Feb 18 peter 12
3723 05 Feb 18 peter 13   The yat library is distributed in the hope that it will be useful,
3723 05 Feb 18 peter 14   but WITHOUT ANY WARRANTY; without even the implied warranty of
3723 05 Feb 18 peter 15   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
3723 05 Feb 18 peter 16   General Public License for more details.
3723 05 Feb 18 peter 17
3723 05 Feb 18 peter 18   You should have received a copy of the GNU General Public License
3723 05 Feb 18 peter 19   along with yat. If not, see <http://www.gnu.org/licenses/>.
3723 05 Feb 18 peter 20 */
3723 05 Feb 18 peter 21
3723 05 Feb 18 peter 22 #include <config.h>
3723 05 Feb 18 peter 23
3723 05 Feb 18 peter 24 #include "Suite.h"
3723 05 Feb 18 peter 25
3723 05 Feb 18 peter 26 #include "yat/statistics/Chi2.h"
3723 05 Feb 18 peter 27
3723 05 Feb 18 peter 28 #include "yat/utility/Matrix.h"
3723 05 Feb 18 peter 29 #include "yat/utility/VectorBase.h"
3723 05 Feb 18 peter 30
3726 09 Feb 18 peter 31 #include <numeric>
3723 05 Feb 18 peter 32
3723 05 Feb 18 peter 33 using namespace theplu::yat;
3723 05 Feb 18 peter 34
3725 08 Feb 18 peter 35 void test1(test::Suite& suite);
3723 05 Feb 18 peter 36
3723 05 Feb 18 peter 37 int main(int argc, char* argv[])
3723 05 Feb 18 peter 38 {
3723 05 Feb 18 peter 39   test::Suite suite(argc, argv);
3723 05 Feb 18 peter 40
3725 08 Feb 18 peter 41   test1(suite);
3723 05 Feb 18 peter 42   return suite.return_value();
3723 05 Feb 18 peter 43 }
3723 05 Feb 18 peter 44
3723 05 Feb 18 peter 45
3725 08 Feb 18 peter 46 void test1(test::Suite& suite)
3723 05 Feb 18 peter 47 {
3723 05 Feb 18 peter 48   suite.out() << "test expected\n";
3723 05 Feb 18 peter 49   statistics::Chi2 chi2;
3723 05 Feb 18 peter 50   utility::Matrix x(2, 3);
3725 08 Feb 18 peter 51   x(0, 0) = 20;
3725 08 Feb 18 peter 52   x(1, 0) = 70;
3725 08 Feb 18 peter 53   x(0, 1) = 50;
3725 08 Feb 18 peter 54   x(1, 1) = 100;
3725 08 Feb 18 peter 55   x(0, 2) = 100;
3725 08 Feb 18 peter 56   x(1, 2) = 200;
3723 05 Feb 18 peter 57
3723 05 Feb 18 peter 58   chi2(x);
3723 05 Feb 18 peter 59   utility::Matrix correct(2, 3);
3723 05 Feb 18 peter 60   for (size_t i=0; i<x.rows(); ++i)
3723 05 Feb 18 peter 61     for (size_t j=0; j<x.columns(); ++j)
3723 05 Feb 18 peter 62       correct(i, j) =
3723 05 Feb 18 peter 63         sum(x.row_const_view(i)) * sum(x.column_const_view(j)) /
3723 05 Feb 18 peter 64         std::accumulate(x.begin(), x.end(), 0.0);
3723 05 Feb 18 peter 65
3723 05 Feb 18 peter 66   if (!suite.equal_matrix(correct, chi2.expected())) {
3723 05 Feb 18 peter 67     suite.err() << "error: incorrect Chi2::expected\n";
3723 05 Feb 18 peter 68     suite.add(false);
3723 05 Feb 18 peter 69   }
3725 08 Feb 18 peter 70
3725 08 Feb 18 peter 71   double p = chi2.p_value();
3725 08 Feb 18 peter 72   if (!suite.equal_fix(p, 0.11692, 0.000001)) {
3725 08 Feb 18 peter 73     suite.err() << "p: " << p << "\n";
3725 08 Feb 18 peter 74     suite.add(false);
3725 08 Feb 18 peter 75   }
3725 08 Feb 18 peter 76
3725 08 Feb 18 peter 77   double mc_p = chi2.p_value(10000);
3725 08 Feb 18 peter 78   if (!suite.equal_fix(mc_p, p, 0.001)) {
3725 08 Feb 18 peter 79     suite.err() << "MC p: " << mc_p << "\n";
3725 08 Feb 18 peter 80     suite.add(false);
3725 08 Feb 18 peter 81   }
3723 05 Feb 18 peter 82 }