test/kendall.cc

Code
Comments
Other
Rev Date Author Line
2649 14 Nov 11 peter 1 // $Id$
2649 14 Nov 11 peter 2
2649 14 Nov 11 peter 3 /*
4089 07 Sep 21 peter 4   Copyright (C) 2011, 2012, 2020, 2021 Peter Johansson
2649 14 Nov 11 peter 5
2649 14 Nov 11 peter 6   This file is part of the yat library, http://dev.thep.lu.se/yat
2649 14 Nov 11 peter 7
2649 14 Nov 11 peter 8   The yat library is free software; you can redistribute it and/or
2649 14 Nov 11 peter 9   modify it under the terms of the GNU General Public License as
2649 14 Nov 11 peter 10   published by the Free Software Foundation; either version 3 of the
2649 14 Nov 11 peter 11   License, or (at your option) any later version.
2649 14 Nov 11 peter 12
2649 14 Nov 11 peter 13   The yat library is distributed in the hope that it will be useful,
2649 14 Nov 11 peter 14   but WITHOUT ANY WARRANTY; without even the implied warranty of
2649 14 Nov 11 peter 15   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
2649 14 Nov 11 peter 16   General Public License for more details.
2649 14 Nov 11 peter 17
2649 14 Nov 11 peter 18   You should have received a copy of the GNU General Public License
2649 14 Nov 11 peter 19   along with yat. If not, see <http://www.gnu.org/licenses/>.
2649 14 Nov 11 peter 20 */
2649 14 Nov 11 peter 21
2881 18 Nov 12 peter 22 #include <config.h>
2881 18 Nov 12 peter 23
2649 14 Nov 11 peter 24 #include "Suite.h"
2649 14 Nov 11 peter 25
3987 26 Aug 20 peter 26 #include "yat/random/random.h"
2649 14 Nov 11 peter 27 #include "yat/statistics/Kendall.h"
3987 26 Aug 20 peter 28 #include "yat/utility/Matrix.h"
4083 27 Aug 21 peter 29 #include "yat/utility/utility.h"
2649 14 Nov 11 peter 30
2758 08 Jul 12 peter 31 #include <gsl/gsl_cdf.h>
2758 08 Jul 12 peter 32
2758 08 Jul 12 peter 33 #include <cmath>
3987 26 Aug 20 peter 34 #include <fstream>
2758 08 Jul 12 peter 35
2649 14 Nov 11 peter 36 using namespace theplu::yat;
2649 14 Nov 11 peter 37 using statistics::Kendall;
2649 14 Nov 11 peter 38
2649 14 Nov 11 peter 39 void test_add(test::Suite&);
2760 08 Jul 12 peter 40 void test_copy(test::Suite&);
2650 14 Nov 11 peter 41 void test_p(test::Suite&);
2650 14 Nov 11 peter 42 void test_p_exact(test::Suite&);
3987 26 Aug 20 peter 43 void test_p_with_ties(test::Suite&);
2650 14 Nov 11 peter 44 void test_score(test::Suite&);
2650 14 Nov 11 peter 45 void test_score_with_ties(test::Suite&);
2649 14 Nov 11 peter 46
2649 14 Nov 11 peter 47 int main(int argc, char* argv[])
2649 14 Nov 11 peter 48 {
2649 14 Nov 11 peter 49   test::Suite suite(argc, argv);
4083 27 Aug 21 peter 50   using std::throw_with_nested;
4083 27 Aug 21 peter 51   try {
4083 27 Aug 21 peter 52     try {
4083 27 Aug 21 peter 53       test_add(suite);
4083 27 Aug 21 peter 54     }
4083 27 Aug 21 peter 55     catch (std::exception& e) {
4083 27 Aug 21 peter 56       throw_with_nested(std::runtime_error("test_add() failed"));
4083 27 Aug 21 peter 57     }
4083 27 Aug 21 peter 58     test_copy(suite);
4083 27 Aug 21 peter 59     test_p(suite);
4110 27 Sep 21 peter 60     test_p_exact(suite);
4110 27 Sep 21 peter 61     test_p_with_ties(suite);
2649 14 Nov 11 peter 62
4083 27 Aug 21 peter 63     try {
4110 27 Sep 21 peter 64       test_score(suite);
4083 27 Aug 21 peter 65     }
4083 27 Aug 21 peter 66     catch (std::exception& e) {
4083 27 Aug 21 peter 67       throw_with_nested(std::runtime_error("test_score() failed"));
4083 27 Aug 21 peter 68     }
4083 27 Aug 21 peter 69
4083 27 Aug 21 peter 70     try {
4110 27 Sep 21 peter 71       test_score_with_ties(suite);
4083 27 Aug 21 peter 72     }
4083 27 Aug 21 peter 73     catch (std::exception& e) {
4083 27 Aug 21 peter 74       throw_with_nested(std::runtime_error("test_score_with_ties() failed"));
4083 27 Aug 21 peter 75     }
4083 27 Aug 21 peter 76   }
4083 27 Aug 21 peter 77   catch (std::exception& e) {
4083 27 Aug 21 peter 78     suite.add(false);
4083 27 Aug 21 peter 79     suite.err() << "error: exception caught: what(): ";
4083 27 Aug 21 peter 80     utility::print_what(e, suite.err());
4083 27 Aug 21 peter 81   }
2649 14 Nov 11 peter 82   return suite.return_value();
2649 14 Nov 11 peter 83 }
2649 14 Nov 11 peter 84
2649 14 Nov 11 peter 85
2649 14 Nov 11 peter 86 void test_add(test::Suite& suite)
2649 14 Nov 11 peter 87 {
2649 14 Nov 11 peter 88   suite.out() << YAT_TEST_PROLOGUE;
4083 27 Aug 21 peter 89   suite.err() << "constructor\n";
2649 14 Nov 11 peter 90   Kendall kendall;
4083 27 Aug 21 peter 91   suite.err() << "::add(0,0)\n";
2649 14 Nov 11 peter 92   kendall.add(0,0);
4083 27 Aug 21 peter 93   suite.err() << "::add(1,1)\n";
2649 14 Nov 11 peter 94   kendall.add(1,1);
4083 27 Aug 21 peter 95   suite.err() << "::score()\n";
4083 27 Aug 21 peter 96   try {
4083 27 Aug 21 peter 97     double score = kendall.score();
4083 27 Aug 21 peter 98     if (!suite.add(suite.equal(score, 1.0)))
4083 27 Aug 21 peter 99       suite.err() << "score not 1.0\n";
4083 27 Aug 21 peter 100   }
4083 27 Aug 21 peter 101   catch (std::exception& e) {
4083 27 Aug 21 peter 102     std::throw_with_nested(std::runtime_error("::score(void) failed"));
4083 27 Aug 21 peter 103   }
2649 14 Nov 11 peter 104 }
2650 14 Nov 11 peter 105
2650 14 Nov 11 peter 106
2760 08 Jul 12 peter 107 void test_copy(test::Suite& suite)
2760 08 Jul 12 peter 108 {
2760 08 Jul 12 peter 109   Kendall k1;
2760 08 Jul 12 peter 110   Kendall k2(k1);
2760 08 Jul 12 peter 111   k1 = k2;
2760 08 Jul 12 peter 112   k2 = k2;
2760 08 Jul 12 peter 113 }
2760 08 Jul 12 peter 114
2760 08 Jul 12 peter 115
2650 14 Nov 11 peter 116 void test_p(test::Suite& suite)
2650 14 Nov 11 peter 117 {
2650 14 Nov 11 peter 118   suite.out() << YAT_TEST_PROLOGUE;
2650 14 Nov 11 peter 119   Kendall kendall;
2650 14 Nov 11 peter 120   // Table 8.1 in Hollander & Wolfe
2650 14 Nov 11 peter 121   kendall.add(44.4, 2.6);
2650 14 Nov 11 peter 122   kendall.add(45.9, 3.1);
2650 14 Nov 11 peter 123   kendall.add(41.9, 2.5);
2650 14 Nov 11 peter 124   kendall.add(53.3, 5.0);
2650 14 Nov 11 peter 125   kendall.add(44.7, 3.6);
2650 14 Nov 11 peter 126   kendall.add(44.1, 4.0);
2650 14 Nov 11 peter 127   kendall.add(50.7, 5.2);
2650 14 Nov 11 peter 128   kendall.add(45.2, 2.8);
2650 14 Nov 11 peter 129   kendall.add(60.1, 3.8);
2650 14 Nov 11 peter 130
2650 14 Nov 11 peter 131   double tau = 16.0 / (9*8/2);
2650 14 Nov 11 peter 132   suite.add(suite.equal(tau, kendall.score()));
2662 21 Nov 11 peter 133   suite.add(suite.equal_fix(0.06, kendall.p_right(true),0.005));
2662 21 Nov 11 peter 134   suite.add(suite.equal_fix(0.04764, kendall.p_right(false),0.00005));
2650 14 Nov 11 peter 135 }
2650 14 Nov 11 peter 136
2650 14 Nov 11 peter 137
2650 14 Nov 11 peter 138 void test_p_exact(test::Suite& suite)
2650 14 Nov 11 peter 139 {
2650 14 Nov 11 peter 140   suite.out() << YAT_TEST_PROLOGUE;
2650 14 Nov 11 peter 141   Kendall kendall;
2650 14 Nov 11 peter 142   for (size_t i=1; i<5; ++i)
2650 14 Nov 11 peter 143     kendall.add(i, i);
2662 21 Nov 11 peter 144   suite.add(suite.equal_fix(0.042, kendall.p_right(true), 0.0005));
2650 14 Nov 11 peter 145   kendall.add(5,5);
2662 21 Nov 11 peter 146   suite.add(suite.equal_fix(0.008, kendall.p_right(true), 0.0005));
2650 14 Nov 11 peter 147   kendall.add(6,6);
2650 14 Nov 11 peter 148   kendall.add(7,7);
2650 14 Nov 11 peter 149   kendall.add(8,1.5);
2662 21 Nov 11 peter 150   suite.add(suite.equal_fix(0.031, kendall.p_right(true), 0.0005));
2650 14 Nov 11 peter 151 }
2650 14 Nov 11 peter 152
2650 14 Nov 11 peter 153
3987 26 Aug 20 peter 154 void test_p_with_ties(test::Suite& suite)
3987 26 Aug 20 peter 155 {
3987 26 Aug 20 peter 156   suite.out() << YAT_TEST_PROLOGUE;
3987 26 Aug 20 peter 157   std::string fn = test::filename("data/kendall.txt");
3987 26 Aug 20 peter 158   std::ifstream is(fn.c_str());
3987 26 Aug 20 peter 159   if (!is.good()) {
3987 26 Aug 20 peter 160     suite.add(false);
3987 26 Aug 20 peter 161     suite.err() << "failed open '" << fn << "'\n";
3987 26 Aug 20 peter 162     return;
3987 26 Aug 20 peter 163   }
3987 26 Aug 20 peter 164   utility::Matrix x(is);
3987 26 Aug 20 peter 165   is.close();
3987 26 Aug 20 peter 166   Kendall kendall;
3987 26 Aug 20 peter 167   for (size_t i=0; i<x.rows(); ++i)
3987 26 Aug 20 peter 168     kendall.add(x(i,0), x(i, 1));
3987 26 Aug 20 peter 169   double score = kendall.score();
3987 26 Aug 20 peter 170   suite.out() << "score: " << score << "\n";
3987 26 Aug 20 peter 171   double p = kendall.p_value();
3987 26 Aug 20 peter 172   suite.out() << "p: " << p << "\n";
4110 27 Sep 21 peter 173   if (p < 0.00001 || p > 0.0001) {
4110 27 Sep 21 peter 174     suite.err() << "error: p not within [0.00001, 0.0001]\n";
3987 26 Aug 20 peter 175     suite.add(false);
3987 26 Aug 20 peter 176   }
3987 26 Aug 20 peter 177 }
3987 26 Aug 20 peter 178
3987 26 Aug 20 peter 179
2650 14 Nov 11 peter 180 void test_score(test::Suite& suite)
2650 14 Nov 11 peter 181 {
2650 14 Nov 11 peter 182   suite.out() << YAT_TEST_PROLOGUE;
2650 14 Nov 11 peter 183   Kendall kendall;
2650 14 Nov 11 peter 184   kendall.add(0,0);
2650 14 Nov 11 peter 185   kendall.add(1,1);
2650 14 Nov 11 peter 186   kendall.add(2,2);
2650 14 Nov 11 peter 187   suite.equal(kendall.score(), 1.0);
2650 14 Nov 11 peter 188   kendall.reset();
2650 14 Nov 11 peter 189   kendall.add(0,2);
2650 14 Nov 11 peter 190   kendall.add(1,1);
2650 14 Nov 11 peter 191   kendall.add(2,0);
2650 14 Nov 11 peter 192   suite.equal(kendall.score(), -1.0);
2650 14 Nov 11 peter 193   kendall.reset();
2650 14 Nov 11 peter 194   kendall.add(0,5);
2650 14 Nov 11 peter 195   kendall.add(1,1);
2650 14 Nov 11 peter 196   kendall.add(2,2);
2650 14 Nov 11 peter 197   kendall.add(3,3);
2650 14 Nov 11 peter 198   suite.equal(kendall.score(), 0.0);
2650 14 Nov 11 peter 199 }
2650 14 Nov 11 peter 200
2650 14 Nov 11 peter 201
2650 14 Nov 11 peter 202 void test_score_with_ties(test::Suite& suite)
2650 14 Nov 11 peter 203 {
2650 14 Nov 11 peter 204   suite.out() << YAT_TEST_PROLOGUE;
2650 14 Nov 11 peter 205   Kendall kendall;
2757 08 Jul 12 peter 206   kendall.add(0,1);
2757 08 Jul 12 peter 207   kendall.add(0,1);
2757 08 Jul 12 peter 208   kendall.add(0,1);
2757 08 Jul 12 peter 209   kendall.add(1,2);
2757 08 Jul 12 peter 210   kendall.add(1,2);
2758 08 Jul 12 peter 211   kendall.add(1,3);
2758 08 Jul 12 peter 212   int nc = 9;
2758 08 Jul 12 peter 213   int nd = 0;
2758 08 Jul 12 peter 214   double correct = (nc-nd) / (std::sqrt(9.0)*std::sqrt(11.0));
2759 08 Jul 12 peter 215   double score = kendall.score();
2759 08 Jul 12 peter 216   suite.out() << "score: " << score << "\n";
2761 08 Jul 12 peter 217   suite.add(suite.equal(correct, score));
2758 08 Jul 12 peter 218   double p = kendall.p_right(false);
2759 08 Jul 12 peter 219   suite.out() << "approx p-value: " << p << "\n";
2758 08 Jul 12 peter 220   int n = kendall.n();
2758 08 Jul 12 peter 221   double v0 = n*(n-1)*(2*n+5);
2758 08 Jul 12 peter 222   double vt = 3*2*11 + 3*2*11;
2761 08 Jul 12 peter 223   double vu = 3*2*11 + 2*1*9;
2761 08 Jul 12 peter 224   double v1 = (6+6)*(6+2)/(2.0*n*(n-1));
3987 26 Aug 20 peter 225   double v2 = (3*2*1+3*2*1)*(3*2*1) / (9.0*n*(n-1)*(n-2));
2758 08 Jul 12 peter 226   double v = (v0-vt-vu)/18 + v1 + v2;
2758 08 Jul 12 peter 227   double z = (nc - nd)/std::sqrt(v);
2758 08 Jul 12 peter 228   correct = gsl_cdf_ugaussian_Q(z);
2761 08 Jul 12 peter 229   suite.add(suite.equal(correct, p));
2759 08 Jul 12 peter 230   p = kendall.p_right(true);
2759 08 Jul 12 peter 231   suite.out() << "exact p-value: " << p << "\n";
2759 08 Jul 12 peter 232   suite.add(suite.equal(36.0/(6*5*4*3*2),p));
2650 14 Nov 11 peter 233 }

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