yat/statistics/Fisher.cc

Code
Comments
Other
Rev Date Author Line
186 07 Oct 04 peter 1 // $Id$
186 07 Oct 04 peter 2
675 10 Oct 06 jari 3 /*
2119 12 Dec 09 peter 4   Copyright (C) 2004 Jari Häkkinen, Peter Johansson
831 27 Mar 07 peter 5   Copyright (C) 2005 Peter Johansson
4359 23 Aug 23 peter 6   Copyright (C) 2006 Jari Häkkinen, Peter Johansson
4359 23 Aug 23 peter 7   Copyright (C) 2007 Peter Johansson
4359 23 Aug 23 peter 8   Copyright (C) 2008 Jari Häkkinen, Peter Johansson
4120 23 Nov 21 peter 9   Copyright (C) 2009, 2010, 2011, 2012, 2013, 2014, 2015, 2018, 2021 Peter Johansson
186 07 Oct 04 peter 10
1437 25 Aug 08 peter 11   This file is part of the yat library, http://dev.thep.lu.se/yat
675 10 Oct 06 jari 12
675 10 Oct 06 jari 13   The yat library is free software; you can redistribute it and/or
675 10 Oct 06 jari 14   modify it under the terms of the GNU General Public License as
1486 09 Sep 08 jari 15   published by the Free Software Foundation; either version 3 of the
675 10 Oct 06 jari 16   License, or (at your option) any later version.
675 10 Oct 06 jari 17
675 10 Oct 06 jari 18   The yat library is distributed in the hope that it will be useful,
675 10 Oct 06 jari 19   but WITHOUT ANY WARRANTY; without even the implied warranty of
675 10 Oct 06 jari 20   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
675 10 Oct 06 jari 21   General Public License for more details.
675 10 Oct 06 jari 22
675 10 Oct 06 jari 23   You should have received a copy of the GNU General Public License
1487 10 Sep 08 jari 24   along with yat. If not, see <http://www.gnu.org/licenses/>.
675 10 Oct 06 jari 25 */
675 10 Oct 06 jari 26
2881 18 Nov 12 peter 27 #include <config.h>
2881 18 Nov 12 peter 28
680 11 Oct 06 jari 29 #include "Fisher.h"
680 11 Oct 06 jari 30 #include "utility.h"
675 10 Oct 06 jari 31
2210 05 Mar 10 peter 32 #include "yat/utility/Exception.h"
2210 05 Mar 10 peter 33
449 15 Dec 05 peter 34 #include <gsl/gsl_cdf.h>
449 15 Dec 05 peter 35 #include <gsl/gsl_randist.h>
449 15 Dec 05 peter 36
1746 23 Jan 09 peter 37 #include <algorithm>
4120 23 Nov 21 peter 38 #include <sstream>
2382 21 Dec 10 peter 39 #include <cassert>
1746 23 Jan 09 peter 40
186 07 Oct 04 peter 41 namespace theplu {
680 11 Oct 06 jari 42 namespace yat {
295 29 Apr 05 peter 43 namespace statistics {
186 07 Oct 04 peter 44
3004 24 Mar 13 peter 45
3298 08 Aug 14 peter 46   Fisher::Fisher(bool yates)
3298 08 Aug 14 peter 47     : a_(0), b_(0), c_(0), d_(0), minimum_size_(10), oddsratio_(1.0),
3298 08 Aug 14 peter 48       yates_(yates)
186 07 Oct 04 peter 49   {
186 07 Oct 04 peter 50   }
186 07 Oct 04 peter 51
186 07 Oct 04 peter 52
777 04 Mar 07 peter 53   Fisher::~Fisher()
777 04 Mar 07 peter 54   {
777 04 Mar 07 peter 55   }
777 04 Mar 07 peter 56
777 04 Mar 07 peter 57
777 04 Mar 07 peter 58   bool Fisher::calculate_p_exact(void) const
777 04 Mar 07 peter 59   {
3004 24 Mar 13 peter 60     return ( a_<minimum_size_ || b_<minimum_size_ ||
3004 24 Mar 13 peter 61              c_<minimum_size_ || d_<minimum_size_);
777 04 Mar 07 peter 62   }
777 04 Mar 07 peter 63
3004 24 Mar 13 peter 64   double Fisher::Chi2(void) const
447 15 Dec 05 peter 65   {
449 15 Dec 05 peter 66     double a,b,c,d;
449 15 Dec 05 peter 67     expected(a,b,c,d);
3298 08 Aug 14 peter 68     if (yates_)
3298 08 Aug 14 peter 69       return yates(a_, a) + yates(b_, b) + yates(c_, c) + yates(d_, d);
3298 08 Aug 14 peter 70
3004 24 Mar 13 peter 71     return (a-a_)*(a-a_)/a + (b-b_)*(b-b_)/b +
449 15 Dec 05 peter 72       (c-c_)*(c-c_)/c + (d-d_)*(d-d_)/d;
449 15 Dec 05 peter 73   }
449 15 Dec 05 peter 74
449 15 Dec 05 peter 75
449 15 Dec 05 peter 76   void Fisher::expected(double& a, double& b, double& c, double& d) const
449 15 Dec 05 peter 77   {
2382 21 Dec 10 peter 78     // use floting point arithmetic to avoid overflow
2382 21 Dec 10 peter 79     double a1 = a_;
2382 21 Dec 10 peter 80     double b1 = b_;
2382 21 Dec 10 peter 81     double c1 = c_;
2382 21 Dec 10 peter 82     double d1 = d_;
2382 21 Dec 10 peter 83     double N = a1 + b1 + c1 + d1;
2382 21 Dec 10 peter 84     a =((a1+b1)*(a1+c1)) / N;
2382 21 Dec 10 peter 85     b =((a1+b1)*(b1+d1)) / N;
2382 21 Dec 10 peter 86     c =((c1+d1)*(a1+c1)) / N;
2382 21 Dec 10 peter 87     d =((c1+d1)*(b1+d1)) / N;
447 15 Dec 05 peter 88   }
447 15 Dec 05 peter 89
718 26 Dec 06 jari 90
3441 23 Nov 15 peter 91   /*
3743 12 Jul 18 peter 92    (n)       n!
3743 12 Jul 18 peter 93    ( )  = --------
3743 12 Jul 18 peter 94    (k)    k!(n-k)!
3743 12 Jul 18 peter 95
3743 12 Jul 18 peter 96
3743 12 Jul 18 peter 97
3441 23 Nov 15 peter 98     ( n )
3743 12 Jul 18 peter 99     (k+1)   k!(n-k)!
3743 12 Jul 18 peter 100     ----- = ------------ = (n-k) / (k+1)
3743 12 Jul 18 peter 101      (n)    (k+1)!(n-k-1)!
3441 23 Nov 15 peter 102      (k)
3441 23 Nov 15 peter 103    */
3441 23 Nov 15 peter 104   double Fisher::choose_ratio(unsigned int n, unsigned int k) const
3441 23 Nov 15 peter 105   {
3743 12 Jul 18 peter 106     assert(k+1 <= n);
3441 23 Nov 15 peter 107     return static_cast<double>(n-k) / (k+1);
3441 23 Nov 15 peter 108   }
3441 23 Nov 15 peter 109
3441 23 Nov 15 peter 110
3441 23 Nov 15 peter 111   double Fisher::hypergeometric_ratio(unsigned int k, unsigned int n1,
3441 23 Nov 15 peter 112                                       unsigned int n2, unsigned int t) const
3441 23 Nov 15 peter 113   {
3743 12 Jul 18 peter 114     // P(X = k+1) / P(X = k) =
3743 12 Jul 18 peter 115     // (choose(n1, k+1) * choose(n2, t-k-1) / choose(n1+n2, t) ) /
3743 12 Jul 18 peter 116     // (choose(n1, k) * choose(n2, t-k) / choose(n1+n2, t) ) =
3743 12 Jul 18 peter 117     // choose(n1, k+1) / choose(n1, k) /( choose(n2, t-k)/choose(n2,t-k-1) ) =
3743 12 Jul 18 peter 118     // choose(n1,k+1)/choose(n1,k) / (choose(n2,(t-k-1)+1)/choose(n2,t-k-1))
3743 12 Jul 18 peter 119     // choose_ratio(n1,k) / choose_ratio(n2, t-k-1)
3743 12 Jul 18 peter 120     // where choose_ratio(a, b) = choose(a, b+1) / choose(a,b)
3441 23 Nov 15 peter 121     return choose_ratio(n1, k) / choose_ratio(n2, t-k-1);
3441 23 Nov 15 peter 122   }
3441 23 Nov 15 peter 123
3441 23 Nov 15 peter 124   double Fisher::lower_tail(unsigned int k, unsigned int n1, unsigned int n2,
3441 23 Nov 15 peter 125                             unsigned int t) const
3441 23 Nov 15 peter 126   {
3441 23 Nov 15 peter 127     double sum = 0;
3441 23 Nov 15 peter 128
3441 23 Nov 15 peter 129     // P(X=k)
3441 23 Nov 15 peter 130     double p0 = gsl_ran_hypergeometric_pdf(k, n1, n2, t);
3441 23 Nov 15 peter 131
3441 23 Nov 15 peter 132     // minimum possible outcome for X is max(0, t-n2)
3441 23 Nov 15 peter 133     unsigned int i = std::max(n2, t)-n2;
3441 23 Nov 15 peter 134     double p = gsl_ran_hypergeometric_pdf(i, n1, n2, t);
3441 23 Nov 15 peter 135     // avoid double dipping P(k)
3441 23 Nov 15 peter 136     for ( ; i<k; ++i) {
3441 23 Nov 15 peter 137       if (p<=p0)
3441 23 Nov 15 peter 138         sum += p;
3441 23 Nov 15 peter 139       else
3441 23 Nov 15 peter 140         break;
3441 23 Nov 15 peter 141
3441 23 Nov 15 peter 142       // calculate p(i+1) using recursive function
3441 23 Nov 15 peter 143       p *= hypergeometric_ratio(i, n1, n2, t);
3441 23 Nov 15 peter 144     }
3441 23 Nov 15 peter 145
3441 23 Nov 15 peter 146     return sum;
3441 23 Nov 15 peter 147   }
3441 23 Nov 15 peter 148
3441 23 Nov 15 peter 149
1271 09 Apr 08 peter 150   unsigned int& Fisher::minimum_size(void)
718 26 Dec 06 jari 151   {
718 26 Dec 06 jari 152     return minimum_size_;
718 26 Dec 06 jari 153   }
718 26 Dec 06 jari 154
718 26 Dec 06 jari 155
1271 09 Apr 08 peter 156   const unsigned int& Fisher::minimum_size(void) const
186 07 Oct 04 peter 157   {
777 04 Mar 07 peter 158     return minimum_size_;
777 04 Mar 07 peter 159   }
777 04 Mar 07 peter 160
777 04 Mar 07 peter 161
1271 09 Apr 08 peter 162   double Fisher::oddsratio(const unsigned int a,
1271 09 Apr 08 peter 163                            const unsigned int b,
1271 09 Apr 08 peter 164                            const unsigned int c,
3004 24 Mar 13 peter 165                            const unsigned int d)
777 04 Mar 07 peter 166   {
186 07 Oct 04 peter 167     // If a column sum or a row sum is zero, the table is nonsense
295 29 Apr 05 peter 168     if ((a==0 || d==0) && (c==0 || b==0)){
4120 23 Nov 21 peter 169       std::ostringstream ss;
4120 23 Nov 21 peter 170       ss << "Table in Fisher is not valid: "
4120 23 Nov 21 peter 171          << a << ", " << b << ", " << c << ", " << d;
4120 23 Nov 21 peter 172       throw utility::runtime_error(ss.str());
186 07 Oct 04 peter 173     }
1746 23 Jan 09 peter 174     a_ = a;
1746 23 Jan 09 peter 175     b_ = b;
1746 23 Jan 09 peter 176     c_ = c;
1746 23 Jan 09 peter 177     d_ = d;
186 07 Oct 04 peter 178
4121 23 Nov 21 peter 179     oddsratio_= (static_cast<double>(a) * static_cast<double>(d)) /
4121 23 Nov 21 peter 180       (static_cast<double>(b) * static_cast<double>(c) );
449 15 Dec 05 peter 181     return oddsratio_;
186 07 Oct 04 peter 182   }
186 07 Oct 04 peter 183
186 07 Oct 04 peter 184
2523 18 Jul 11 peter 185   double Fisher::oddsratio(void) const
2523 18 Jul 11 peter 186   {
2523 18 Jul 11 peter 187     return oddsratio_;
2523 18 Jul 11 peter 188   }
2523 18 Jul 11 peter 189
2523 18 Jul 11 peter 190
3004 24 Mar 13 peter 191   double Fisher::p_left(void) const
186 07 Oct 04 peter 192   {
3004 24 Mar 13 peter 193     if (!calculate_p_exact()) {
3004 24 Mar 13 peter 194       if (oddsratio_>1.0)
3004 24 Mar 13 peter 195         return 1.0-p_value_approximative()/2.0;
3004 24 Mar 13 peter 196       return p_value_approximative()/2.0;
3004 24 Mar 13 peter 197     }
3441 23 Nov 15 peter 198     return p_left_exact();
3441 23 Nov 15 peter 199   }
3441 23 Nov 15 peter 200
3441 23 Nov 15 peter 201
3441 23 Nov 15 peter 202   double Fisher::p_left_exact(void) const
3441 23 Nov 15 peter 203   {
3004 24 Mar 13 peter 204     // check for overflow
3004 24 Mar 13 peter 205     assert(c_ <= c_+d_ && d_ <= c_+d_);
3004 24 Mar 13 peter 206     assert(a_ <= a_+b_ && b_ <= a_+b_);
3004 24 Mar 13 peter 207     assert(a_ <= a_+c_ && c_ <= a_+c_);
3004 24 Mar 13 peter 208
3004 24 Mar 13 peter 209     return cdf_hypergeometric_P(a_, a_+b_, c_+d_, a_+c_);
3004 24 Mar 13 peter 210   }
3004 24 Mar 13 peter 211
3004 24 Mar 13 peter 212
3004 24 Mar 13 peter 213   double Fisher::p_value(void) const
3004 24 Mar 13 peter 214   {
1746 23 Jan 09 peter 215     if (calculate_p_exact())
1746 23 Jan 09 peter 216       return p_value_exact();
1746 23 Jan 09 peter 217     return p_value_approximative();
777 04 Mar 07 peter 218   }
447 15 Dec 05 peter 219
777 04 Mar 07 peter 220
3004 24 Mar 13 peter 221   double Fisher::p_right(void) const
777 04 Mar 07 peter 222   {
1746 23 Jan 09 peter 223     if (!calculate_p_exact()) {
1746 23 Jan 09 peter 224       if (oddsratio_<1.0)
1746 23 Jan 09 peter 225         return 1.0-p_value_approximative()/2.0;
1746 23 Jan 09 peter 226       return p_value_approximative()/2.0;
1746 23 Jan 09 peter 227     }
3441 23 Nov 15 peter 228     return p_right_exact();
3441 23 Nov 15 peter 229   }
3441 23 Nov 15 peter 230
3441 23 Nov 15 peter 231
3441 23 Nov 15 peter 232   double Fisher::p_right_exact(void) const
3441 23 Nov 15 peter 233   {
2382 21 Dec 10 peter 234     // check for overflow
2382 21 Dec 10 peter 235     assert(c_ <= c_+d_ && d_ <= c_+d_);
2382 21 Dec 10 peter 236     assert(a_ <= a_+b_ && b_ <= a_+b_);
2382 21 Dec 10 peter 237     assert(a_ <= a_+c_ && c_ <= a_+c_);
2382 21 Dec 10 peter 238
3004 24 Mar 13 peter 239     return cdf_hypergeometric_P(c_, c_+d_, a_+b_, a_+c_);
447 15 Dec 05 peter 240   }
447 15 Dec 05 peter 241
447 15 Dec 05 peter 242
3004 24 Mar 13 peter 243   double Fisher::p_value_one_sided(void) const
447 15 Dec 05 peter 244   {
3004 24 Mar 13 peter 245     return p_right();
3004 24 Mar 13 peter 246   }
3004 24 Mar 13 peter 247
3004 24 Mar 13 peter 248
3004 24 Mar 13 peter 249   double Fisher::p_value_approximative(void) const
3004 24 Mar 13 peter 250   {
449 15 Dec 05 peter 251     return gsl_cdf_chisq_Q(Chi2(), 1.0);
447 15 Dec 05 peter 252   }
447 15 Dec 05 peter 253
3004 24 Mar 13 peter 254
3004 24 Mar 13 peter 255   double Fisher::p_value_exact(void) const
3004 24 Mar 13 peter 256   {
3441 23 Nov 15 peter 257     // The hypergeometric distribution is unimodal (only one peak) and the
3743 12 Jul 18 peter 258     // peak (mode) occurs at: floor((t+1)*(n1+1)/(n+1))
449 15 Dec 05 peter 259
3743 12 Jul 18 peter 260     double mode = std::floor((a_+c_+1.0)*(a_+b_+1.0)/(a_+b_+c_+d_+1.0));
3441 23 Nov 15 peter 261     if (a_ >= mode)
3441 23 Nov 15 peter 262       return p_value_exact(a_, a_+b_, c_+d_, a_+c_);
3441 23 Nov 15 peter 263     return p_value_exact(c_, c_+d_, a_+b_, a_+c_);
3441 23 Nov 15 peter 264   }
449 15 Dec 05 peter 265
3441 23 Nov 15 peter 266
3441 23 Nov 15 peter 267   /*
3441 23 Nov 15 peter 268     a | b | n1
3441 23 Nov 15 peter 269     c | d | n2
3441 23 Nov 15 peter 270     -------------------
3441 23 Nov 15 peter 271     t |   | n
3441 23 Nov 15 peter 272
3441 23 Nov 15 peter 273
3441 23 Nov 15 peter 274     The p is calculated as the sum of all P(x) for all x such that
3441 23 Nov 15 peter 275     P(x)<=P(k)
3441 23 Nov 15 peter 276   */
3441 23 Nov 15 peter 277   double Fisher::p_value_exact(unsigned int k, unsigned int n1,
3441 23 Nov 15 peter 278                                unsigned int n2, unsigned int t) const
3441 23 Nov 15 peter 279   {
3441 23 Nov 15 peter 280     assert(k<=t);
3441 23 Nov 15 peter 281     assert(t<=n1+n2);
3441 23 Nov 15 peter 282     assert(n1);
3441 23 Nov 15 peter 283     assert(n2);
3441 23 Nov 15 peter 284
3743 12 Jul 18 peter 285     // special case when k=0 and mode (peak of distribution) is at
3743 12 Jul 18 peter 286     // zero as well. If mode is not zero 2x2 table should have been
3743 12 Jul 18 peter 287     // mirrored before calling this function.
3743 12 Jul 18 peter 288     if (k == 0) {
3743 12 Jul 18 peter 289       assert(std::floor((a_+c_+1.0)*(a_+b_+1.0)/(a_+b_+c_+d_+1.0)) == 0.0);
3743 12 Jul 18 peter 290       return 1.0;
3743 12 Jul 18 peter 291     }
3743 12 Jul 18 peter 292
3441 23 Nov 15 peter 293     assert(k);
3441 23 Nov 15 peter 294     // P(X >= k) = P(X > k-1)
3441 23 Nov 15 peter 295     double sum = gsl_cdf_hypergeometric_Q(k-1, n1, n2, t);
3441 23 Nov 15 peter 296     // calculate the other tail
3441 23 Nov 15 peter 297     sum += lower_tail(k, n1, n2, t);
3441 23 Nov 15 peter 298     return sum;
186 07 Oct 04 peter 299   }
186 07 Oct 04 peter 300
3298 08 Aug 14 peter 301
3298 08 Aug 14 peter 302   double Fisher::yates(unsigned int o, double e) const
3298 08 Aug 14 peter 303   {
3298 08 Aug 14 peter 304     double x = std::abs(o - e) - 0.5;
3298 08 Aug 14 peter 305     return x*x/e;
3298 08 Aug 14 peter 306   }
3298 08 Aug 14 peter 307
683 11 Oct 06 jari 308 }}} // of namespace statistics, yat, and theplu

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