source: trunk/yat/statistics/Fisher.cc @ 3455

Last change on this file since 3455 was 3455, checked in by Peter, 7 years ago

update copyright years

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