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

Last change on this file since 1486 was 1486, checked in by Jari Häkkinen, 13 years ago

Addresses #436.

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 3.7 KB
Line 
1// $Id: Fisher.cc 1486 2008-09-09 21:17:19Z jari $
2
3/*
4  Copyright (C) 2004 Jari Häkkinen, Peter Johansson
5  Copyright (C) 2005 Peter Johansson
6  Copyright (C) 2006, 2007 Jari Häkkinen, Peter Johansson
7  Copyright (C) 2008 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 this program; if not, write to the Free Software
23  Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA
24  02111-1307, USA.
25*/
26
27#include "Fisher.h"
28#include "utility.h"
29
30#include <gsl/gsl_cdf.h>
31#include <gsl/gsl_randist.h>
32
33namespace theplu {
34namespace yat {
35namespace statistics {
36
37 
38  Fisher::Fisher()
39    : a_(0), b_(0), c_(0), d_(0), oddsratio_(1.0)
40  {
41  }
42
43
44  Fisher::~Fisher()
45  {
46  }
47
48
49  bool Fisher::calculate_p_exact(void) const
50  {
51    return ( a_<minimum_size_ || b_<minimum_size_ || 
52             c_<minimum_size_ || d_<minimum_size_); 
53  }
54
55  double Fisher::Chi2() const
56  {
57    double a,b,c,d;
58    expected(a,b,c,d);
59    return (a-a_)*(a-a_)/a + (b-b_)*(b-b_)/b + 
60      (c-c_)*(c-c_)/c + (d-d_)*(d-d_)/d;
61  }
62
63
64  void Fisher::expected(double& a, double& b, double& c, double& d) const
65  {
66    double N = a_+b_+c_+d_;
67    a =((a_+b_)*(a_+c_)) / N;
68    b =((a_+b_)*(b_+d_)) / N;
69    c =((c_+d_)*(a_+c_)) / N;
70    d =((c_+d_)*(b_+d_)) / N;
71  }
72
73
74  unsigned int& Fisher::minimum_size(void)
75  {
76    return minimum_size_;
77  }
78
79
80  const unsigned int& Fisher::minimum_size(void) const
81  {
82    return minimum_size_;
83  }
84
85
86  double Fisher::oddsratio(const unsigned int a,
87                           const unsigned int b,
88                           const unsigned int c,
89                           const unsigned int d) 
90  {
91    // If a column sum or a row sum is zero, the table is nonsense
92    if ((a==0 || d==0) && (c==0 || b==0)){
93      throw std::runtime_error("runtime_error: Table in Fisher is not valid\n");
94    }
95
96    oddsratio_=(a*d)/(b*d);
97    return oddsratio_;
98  }
99
100
101  double Fisher::p_value() const
102  {
103    if (oddsratio_>1.0)
104      return 2*p_value_one_sided();
105    if (oddsratio_<1.0){
106      // If oddsratio is less than unity
107      // Two-sided p-value is
108      // P(X <= oddsratio) + P(X >= 1/oddsratio) =
109      // P(X <= oddsratio) + P(X <= oddsratio)
110      // 2 * (P(X < oddsratio) + P(X = oddsratio))
111      // 2 * (1 - P(X >= oddsratio) + P(X = oddsratio))
112      if (calculate_p_exact())
113        return 2*(1-p_value_one_sided()+
114                  gsl_ran_hypergeometric_pdf(a_,a_+b_,c_+d_, a_+c_)
115                  );
116      // if p is calculated approximatively correction is not needed
117      return 2*(1-p_value_one_sided());
118
119    }
120    return 1.0;
121  }
122
123
124  double Fisher::p_value_one_sided() const
125  {
126    if ( calculate_p_exact() )
127      return p_value_exact();
128    return p_value_approximative();
129  }
130
131
132  double Fisher::p_value_approximative() const
133  {
134    return gsl_cdf_chisq_Q(Chi2(), 1.0);
135  }
136
137  double Fisher::p_value_exact() const
138  { 
139    // Since the calculation is symmetric and cdf_hypergeometric_P
140    // loops up to k we choose the smallest number to be k and mirror
141    // the matrix. This choice makes the p-value two-sided.
142
143    if (a_<b_ && a_<c_ && a_<d_)
144      return statistics::cdf_hypergeometric_P(a_,a_+b_,c_+d_,a_+c_);
145    else if (b_<a_ && b_<c_ && b_<d_)
146      return  statistics::cdf_hypergeometric_P(b_,a_+b_,c_+d_,b_+d_);
147    else if (c_<a_ && c_<b_ && c_<d_)
148      return  statistics::cdf_hypergeometric_P(c_,c_+d_,a_+b_,a_+c_);
149
150    return statistics::cdf_hypergeometric_P(d_,c_+d_,a_+b_,b_+d_);
151
152  }
153
154
155}}} // of namespace statistics, yat, and theplu
Note: See TracBrowser for help on using the repository browser.