test/fisher.cc

Code
Comments
Other
Rev Date Author Line
1621 10 Nov 08 peter 1 // $Id$
1621 10 Nov 08 peter 2
1621 10 Nov 08 peter 3 /*
3743 12 Jul 18 peter 4   Copyright (C) 2008, 2009, 2010, 2012, 2013, 2014, 2015, 2018 Peter Johansson
1621 10 Nov 08 peter 5
1621 10 Nov 08 peter 6   This file is part of the yat library, http://dev.thep.lu.se/yat
1621 10 Nov 08 peter 7
1621 10 Nov 08 peter 8   The yat library is free software; you can redistribute it and/or
1621 10 Nov 08 peter 9   modify it under the terms of the GNU General Public License as
1621 10 Nov 08 peter 10   published by the Free Software Foundation; either version 3 of the
1621 10 Nov 08 peter 11   License, or (at your option) any later version.
1621 10 Nov 08 peter 12
1621 10 Nov 08 peter 13   The yat library is distributed in the hope that it will be useful,
1621 10 Nov 08 peter 14   but WITHOUT ANY WARRANTY; without even the implied warranty of
1621 10 Nov 08 peter 15   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
1621 10 Nov 08 peter 16   General Public License for more details.
1621 10 Nov 08 peter 17
1621 10 Nov 08 peter 18   You should have received a copy of the GNU General Public License
1621 10 Nov 08 peter 19   along with yat. If not, see <http://www.gnu.org/licenses/>.
1621 10 Nov 08 peter 20 */
1621 10 Nov 08 peter 21
2881 18 Nov 12 peter 22 #include <config.h>
2881 18 Nov 12 peter 23
1621 10 Nov 08 peter 24 #include "Suite.h"
1621 10 Nov 08 peter 25
1621 10 Nov 08 peter 26 #include "yat/statistics/Fisher.h"
1621 10 Nov 08 peter 27
3743 12 Jul 18 peter 28 #include <gsl/gsl_cdf.h>
3743 12 Jul 18 peter 29 #include <gsl/gsl_randist.h>
2382 21 Dec 10 peter 30 #include <climits>
2382 21 Dec 10 peter 31
1624 12 Nov 08 peter 32 using namespace theplu::yat;
1624 12 Nov 08 peter 33 void test_p_value(test::Suite&);
1624 12 Nov 08 peter 34 void test_p_value_approximative(test::Suite&);
1624 12 Nov 08 peter 35 void test_p_value_exact(test::Suite&);
2382 21 Dec 10 peter 36 void test_large_numbers(test::Suite&);
3298 08 Aug 14 peter 37 void test_yates(test::Suite&);
3743 12 Jul 18 peter 38 void test_cdf_hypergeometric_Q(test::Suite&);
3743 12 Jul 18 peter 39 void test_ticket908(test::Suite&);
1624 12 Nov 08 peter 40
1621 10 Nov 08 peter 41 int main(int argc, char* argv[])
2706 13 Mar 12 peter 42 {
1621 10 Nov 08 peter 43   test::Suite suite(argc, argv);
1621 10 Nov 08 peter 44
1621 10 Nov 08 peter 45   statistics::Fisher f;
1621 10 Nov 08 peter 46   if (!suite.equal(f.oddsratio(1,4,5,1), 0.05)) {
1621 10 Nov 08 peter 47     suite.add(false);
1621 10 Nov 08 peter 48     suite.err() << "oddsratio failed\n";
1621 10 Nov 08 peter 49   }
1621 10 Nov 08 peter 50   double a, b, c, d;
1621 10 Nov 08 peter 51   f.expected(a,b,c,d);
1621 10 Nov 08 peter 52   if (!suite.equal(a, 30.0/11.0)) {
1621 10 Nov 08 peter 53     suite.add(false);
1621 10 Nov 08 peter 54     suite.err() << "expected a failed\n";
1621 10 Nov 08 peter 55   }
1621 10 Nov 08 peter 56   if (!suite.equal(b, 25.0/11.0)) {
1621 10 Nov 08 peter 57     suite.add(false);
1621 10 Nov 08 peter 58     suite.err() << "expected b failed\n";
1621 10 Nov 08 peter 59   }
1621 10 Nov 08 peter 60   if (!suite.equal(c, 36.0/11.0)) {
1621 10 Nov 08 peter 61     suite.add(false);
1621 10 Nov 08 peter 62     suite.err() << "expected c failed\n";
1621 10 Nov 08 peter 63   }
1621 10 Nov 08 peter 64   if (!suite.equal(d, 30.0/11.0)) {
1621 10 Nov 08 peter 65     suite.add(false);
1621 10 Nov 08 peter 66     suite.err() << "expected d failed\n";
1621 10 Nov 08 peter 67   }
2706 13 Mar 12 peter 68   if (!suite.equal(f.Chi2(), 4.4122222222222222222222)) {
1621 10 Nov 08 peter 69     suite.add(false);
1621 10 Nov 08 peter 70     suite.err() << "Chi2 failed\n";
1621 10 Nov 08 peter 71   }
1621 10 Nov 08 peter 72   if (!suite.equal(f.minimum_size(),10)) {
1621 10 Nov 08 peter 73     suite.add(false);
1621 10 Nov 08 peter 74     suite.err() << "minimum_size failed\n";
1621 10 Nov 08 peter 75   }
1624 12 Nov 08 peter 76   test_p_value(suite);
2382 21 Dec 10 peter 77   test_large_numbers(suite);
3298 08 Aug 14 peter 78   test_yates(suite);
3743 12 Jul 18 peter 79   test_cdf_hypergeometric_Q(suite);
3743 12 Jul 18 peter 80   test_ticket908(suite);
1621 10 Nov 08 peter 81   return suite.return_value();
1621 10 Nov 08 peter 82 }
1624 12 Nov 08 peter 83
1624 12 Nov 08 peter 84
2382 21 Dec 10 peter 85 void test_large_numbers(test::Suite& suite)
2382 21 Dec 10 peter 86 {
2382 21 Dec 10 peter 87   // skip test if unsigned int is 16 bit
2382 21 Dec 10 peter 88   if ((UINT_MAX >> 16) == 0) {
2382 21 Dec 10 peter 89     suite.out() << "skipping test_large_numbers\n";
2382 21 Dec 10 peter 90     return;
2382 21 Dec 10 peter 91   }
2706 13 Mar 12 peter 92
2382 21 Dec 10 peter 93   statistics::Fisher f;
2382 21 Dec 10 peter 94   double oddsratio = f.oddsratio(1166,63326825-1166,1095,66074759-1095);
2382 21 Dec 10 peter 95   if (oddsratio<0.5 || oddsratio>2) {
2382 21 Dec 10 peter 96     suite.err() << "oddsratio: " << oddsratio << "\n";
2382 21 Dec 10 peter 97     suite.err() << "expected ~ 1\n";
2706 13 Mar 12 peter 98     suite.add(false);
2382 21 Dec 10 peter 99   }
2382 21 Dec 10 peter 100   suite.add(suite.equal_fix(f.p_value(), 0.0123, 0.0001));
3004 24 Mar 13 peter 101   f.p_left();
3004 24 Mar 13 peter 102   f.p_right();
2382 21 Dec 10 peter 103 }
2382 21 Dec 10 peter 104
2382 21 Dec 10 peter 105
1624 12 Nov 08 peter 106 void test_p_value(test::Suite& suite)
1624 12 Nov 08 peter 107 {
1624 12 Nov 08 peter 108   test_p_value_exact(suite);
1624 12 Nov 08 peter 109   test_p_value_approximative(suite);
1624 12 Nov 08 peter 110 }
1624 12 Nov 08 peter 111
1624 12 Nov 08 peter 112
1624 12 Nov 08 peter 113 void test_p_value_approximative(test::Suite& suite)
1624 12 Nov 08 peter 114 {
1624 12 Nov 08 peter 115   suite.err() << "testing p_value_approximative\n";
1624 12 Nov 08 peter 116   statistics::Fisher f;
1624 12 Nov 08 peter 117   f.minimum_size() = 0;
1624 12 Nov 08 peter 118   f.oddsratio(10,20,20,50);
3004 24 Mar 13 peter 119   // oddsratio 1.25 > 1 so p = 2*p_right
3004 24 Mar 13 peter 120   suite.add(suite.equal(f.p_value(), 2*f.p_right()));
1624 12 Nov 08 peter 121   f.oddsratio(10,20,10,20);
1624 12 Nov 08 peter 122   suite.add(suite.equal(f.p_value(), 1.0));
3004 24 Mar 13 peter 123   suite.add(suite.equal(f.p_right(), 0.5));
3004 24 Mar 13 peter 124   suite.add(suite.equal(f.p_left(), 0.5));
1624 12 Nov 08 peter 125 }
1624 12 Nov 08 peter 126
1624 12 Nov 08 peter 127
1624 12 Nov 08 peter 128 void test_p_value_exact(test::Suite& suite)
1624 12 Nov 08 peter 129 {
1624 12 Nov 08 peter 130   suite.err() << "test p_value_exact\n";
1624 12 Nov 08 peter 131   statistics::Fisher f;
1624 12 Nov 08 peter 132   f.minimum_size() = 1000;
1624 12 Nov 08 peter 133   f.oddsratio(10,20,10,20);
3441 23 Nov 15 peter 134   suite.add(suite.equal(f.p_value(), 1.0, 20));
1624 12 Nov 08 peter 135
1624 12 Nov 08 peter 136   f.oddsratio(10, 20, 20, 200);
1748 26 Jan 09 peter 137   suite.add(suite.equal_fix(f.p_value(), 0.000811906062767622,1e-16));
3004 24 Mar 13 peter 138   suite.add(suite.equal_fix(f.p_right(), 0.000811906062767622,1e-16));
1624 12 Nov 08 peter 139
1624 12 Nov 08 peter 140   // testing symmetry
1624 12 Nov 08 peter 141   statistics::Fisher f2;
1624 12 Nov 08 peter 142   f2.minimum_size() = 1000;
1624 12 Nov 08 peter 143   f2.oddsratio(20, 200, 10, 20);
1624 12 Nov 08 peter 144   suite.add(suite.equal(f2.p_value(), f.p_value()));
3004 24 Mar 13 peter 145   suite.add(suite.equal(f2.p_left(), f.p_right()));
1624 12 Nov 08 peter 146
1624 12 Nov 08 peter 147   f.oddsratio(1, 1, 1, 2);
1746 23 Jan 09 peter 148   suite.add(suite.equal(f.p_value(), 1.0, 2));
3004 24 Mar 13 peter 149   suite.add(suite.equal(f.p_right(), 0.7));
3004 24 Mar 13 peter 150   suite.add(suite.equal(f.p_left(), 0.9));
1624 12 Nov 08 peter 151
1624 12 Nov 08 peter 152   f.oddsratio(1, 1, 2, 1);
1746 23 Jan 09 peter 153   suite.add(suite.equal(f.p_value(), 1.0, 2));
3004 24 Mar 13 peter 154   suite.add(suite.equal(f.p_right(), 0.9));
3004 24 Mar 13 peter 155   suite.add(suite.equal(f.p_left(), 0.7));
1624 12 Nov 08 peter 156 }
3298 08 Aug 14 peter 157
3298 08 Aug 14 peter 158
3298 08 Aug 14 peter 159 void test_yates(test::Suite& suite)
3298 08 Aug 14 peter 160 {
3298 08 Aug 14 peter 161   statistics::Fisher f;
3298 08 Aug 14 peter 162   statistics::Fisher f2(true);
3298 08 Aug 14 peter 163
3298 08 Aug 14 peter 164   f.oddsratio(5,10,10,10);
3298 08 Aug 14 peter 165   f2.oddsratio(5,10,10,10);
3298 08 Aug 14 peter 166   f.minimum_size() = 0;
3298 08 Aug 14 peter 167   f2.minimum_size() = 0;
3298 08 Aug 14 peter 168   double p = f.p_value();
3298 08 Aug 14 peter 169   double p_yates = f2.p_value();
3298 08 Aug 14 peter 170   if (p==p_yates) {
3298 08 Aug 14 peter 171     suite.add(false);
3298 08 Aug 14 peter 172     suite.err() << "p-value: " << p << "\n";
3298 08 Aug 14 peter 173     suite.err() << "p-value yates's corrected: " << p_yates << "\n";
3298 08 Aug 14 peter 174   }
3298 08 Aug 14 peter 175 }
3743 12 Jul 18 peter 176
3743 12 Jul 18 peter 177
3743 12 Jul 18 peter 178 void test_cdf_hypergeometric_Q(test::Suite& suite)
3743 12 Jul 18 peter 179 {
3743 12 Jul 18 peter 180   suite.out() << "\ntest " << __func__ << "\n";
3743 12 Jul 18 peter 181   unsigned int n1 = 16;
3743 12 Jul 18 peter 182   unsigned int n2 = 3;
3743 12 Jul 18 peter 183   unsigned int t = 16;
3743 12 Jul 18 peter 184   for (unsigned int k=0; k<=t; ++k) {
3743 12 Jul 18 peter 185     // P(X=k)
3743 12 Jul 18 peter 186     double p = gsl_ran_hypergeometric_pdf(k, n1, n2, t);
3743 12 Jul 18 peter 187     // P(X<=k)
3743 12 Jul 18 peter 188     double P = gsl_cdf_hypergeometric_P(k, n1, n2, t);
3743 12 Jul 18 peter 189     // P(X>k)
3743 12 Jul 18 peter 190     double Q = gsl_cdf_hypergeometric_Q(k, n1, n2, t);
3743 12 Jul 18 peter 191     suite.out() << k << " " << p << " " << P << " " << Q << "\n";
3743 12 Jul 18 peter 192     double sum = 0;
3743 12 Jul 18 peter 193     for (unsigned int x=0; x<=k; ++x)
3743 12 Jul 18 peter 194       sum += gsl_ran_hypergeometric_pdf(x, n1, n2, t);
3743 12 Jul 18 peter 195     if (!suite.equal(sum, P, 20)) {
3743 12 Jul 18 peter 196       suite.add(false);
3743 12 Jul 18 peter 197       suite.err() << "error: gsl_cdf_hypergeometric_P for k=" << k << "\n";
3743 12 Jul 18 peter 198     }
3743 12 Jul 18 peter 199     if (!suite.equal(Q, 1-P, 1000)) {
3743 12 Jul 18 peter 200       suite.add(false);
3743 12 Jul 18 peter 201       suite.err() << "failed: Q = 1-P for k=" << k << "\n";
3743 12 Jul 18 peter 202     }
3743 12 Jul 18 peter 203   }
3743 12 Jul 18 peter 204 }
3743 12 Jul 18 peter 205
3743 12 Jul 18 peter 206
3743 12 Jul 18 peter 207 // test bug #908 reported in http://dev.thep.lu.se/yat/ticket/908
3743 12 Jul 18 peter 208 void test_ticket908(test::Suite& suite)
3743 12 Jul 18 peter 209 {
3743 12 Jul 18 peter 210   suite.out() << "===\ntesting ticket 908\n";
3743 12 Jul 18 peter 211   statistics::Fisher f;
3743 12 Jul 18 peter 212   int a1 = 13;
3743 12 Jul 18 peter 213   int a2 = 3;
3743 12 Jul 18 peter 214   int b1 = 3;
3743 12 Jul 18 peter 215   int b2 = 0;
3743 12 Jul 18 peter 216
3743 12 Jul 18 peter 217   double correct_p = 0;
3743 12 Jul 18 peter 218   unsigned int n1 = a1 + a2;
3743 12 Jul 18 peter 219   unsigned int n2 = b1 + b2;
3743 12 Jul 18 peter 220   unsigned int t = a1+b1;
3743 12 Jul 18 peter 221   for (unsigned int k=0; k<=n1; ++k) {
3743 12 Jul 18 peter 222     double P = gsl_ran_hypergeometric_pdf(k, n1, n2, t);
3743 12 Jul 18 peter 223     if (P <=   gsl_ran_hypergeometric_pdf(a1, n1, n2, t)) {
3743 12 Jul 18 peter 224       suite.out() << "counting " << k << " " << P << "\n";
3743 12 Jul 18 peter 225       correct_p += P;
3743 12 Jul 18 peter 226     }
3743 12 Jul 18 peter 227   }
3743 12 Jul 18 peter 228   statistics::Fisher fisher;
3743 12 Jul 18 peter 229   std::vector<double> p;
3743 12 Jul 18 peter 230   fisher.oddsratio(a1, a2, b1, b2);
3743 12 Jul 18 peter 231   p.push_back(fisher.p_value());
3743 12 Jul 18 peter 232   fisher.oddsratio(a2, a1, b2, b1);
3743 12 Jul 18 peter 233   p.push_back(fisher.p_value());
3743 12 Jul 18 peter 234   fisher.oddsratio(b1, b2, a1, a2);
3743 12 Jul 18 peter 235   p.push_back(fisher.p_value());
3743 12 Jul 18 peter 236   fisher.oddsratio(b2, b1, a2, a1);
3743 12 Jul 18 peter 237   p.push_back(fisher.p_value());
3743 12 Jul 18 peter 238   for (size_t i=0; i<p.size(); ++i) {
3743 12 Jul 18 peter 239     if (!suite.add(suite.equal(correct_p, p[i], 20)))
3743 12 Jul 18 peter 240       suite.err() << "error: incorrect p: " << p[i] << " for case "
3743 12 Jul 18 peter 241                   << i << "; expected: " << correct_p << "\n";
3743 12 Jul 18 peter 242   }
3743 12 Jul 18 peter 243 }

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