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

Last change on this file since 760 was 760, checked in by Peter, 15 years ago

throw exception when detecting invalid 2x2 table. Fixes #26

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 4.9 KB
Line 
1// $Id: Fisher.cc 760 2007-02-20 09:58:33Z peter $
2
3/*
4  Copyright (C) The authors contributing to this file.
5
6  This file is part of the yat library, http://lev.thep.lu.se/trac/yat
7
8  The yat library is free software; you can redistribute it and/or
9  modify it under the terms of the GNU General Public License as
10  published by the Free Software Foundation; either version 2 of the
11  License, or (at your option) any later version.
12
13  The yat library is distributed in the hope that it will be useful,
14  but WITHOUT ANY WARRANTY; without even the implied warranty of
15  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
16  General Public License for more details.
17
18  You should have received a copy of the GNU General Public License
19  along with this program; if not, write to the Free Software
20  Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA
21  02111-1307, USA.
22*/
23
24#include "Fisher.h"
25#include "Score.h"
26#include "utility.h"
27#include "yat/classifier/DataLookupWeighted1D.h"
28#include "yat/classifier/Target.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(bool b)
39    : Score(b), a_(0), b_(0), c_(0), d_(0), oddsratio_(1.0), value_cutoff_(0)
40  {
41  }
42
43
44  double Fisher::Chi2() const
45  {
46    double a,b,c,d;
47    expected(a,b,c,d);
48    return (a-a_)*(a-a_)/a + (b-b_)*(b-b_)/b + 
49      (c-c_)*(c-c_)/c + (d-d_)*(d-d_)/d;
50  }
51
52
53  void Fisher::expected(double& a, double& b, double& c, double& d) const
54  {
55    double N = a_+b_+c_+d_;
56    a =((a_+b_)*(a_+c_)) / N;
57    b =((a_+b_)*(b_+d_)) / N;
58    c =((c_+d_)*(a_+c_)) / N;
59    d =((c_+d_)*(b_+d_)) / N;
60  }
61
62
63  u_int& Fisher::minimum_size(void)
64  {
65    return minimum_size_;
66  }
67
68
69  double Fisher::oddsratio(const double a,
70                           const double b,
71                           const double c,
72                           const double d) 
73  {
74    // If a column sum or a row sum is zero, the table is nonsense
75    if ((a==0 || d==0) && (c==0 || b==0)){
76      throw std::runtime_error("runtime_error: Table in Fisher is not valid\n");
77    }
78
79    oddsratio_=(a*d)/(b*d);
80    if (absolute_ && oddsratio_<1)
81      return 1/oddsratio_;
82    return oddsratio_;
83  }
84
85
86  double Fisher::p_value() const
87  {
88    double p=1;
89    if ( ( a_<minimum_size_ || b_<minimum_size_ || 
90           c_<minimum_size_ || d_<minimum_size_) && !weighted_ ){
91
92      p=p_value_exact();
93    }
94    else
95      p=p_value_approximative();
96
97    if (!absolute_){
98      p=p/2;
99      if (oddsratio_<0.5){
100        // last term because >= not equal to !(<=)
101        u_int a = static_cast<u_int>(a_);
102        u_int b = static_cast<u_int>(b_);
103        u_int c = static_cast<u_int>(c_);
104        u_int d = static_cast<u_int>(d_);
105        return 1-p+gsl_ran_hypergeometric_pdf(a, a+b, c+d, a+c);
106      }
107    }
108    return p;
109  }
110
111
112  double Fisher::p_value_approximative() const
113  {
114    return gsl_cdf_chisq_Q(Chi2(), 1.0);
115  }
116
117  double Fisher::p_value_exact() const
118  { 
119    u_int a = static_cast<u_int>(a_);
120    u_int b = static_cast<u_int>(b_);
121    u_int c = static_cast<u_int>(c_);
122    u_int d = static_cast<u_int>(d_);
123   
124    // Since the calculation is symmetric and cdf_hypergeometric_P
125    // loops up to k we choose the smallest number to be k and mirror
126    // the matrix. This choice makes the p-value two-sided.
127
128    if (a<b && a<c && a<d)
129      return statistics::cdf_hypergeometric_P(a,a+b,c+d,a+c);
130    else if (b<a && b<c && b<d)
131      return  statistics::cdf_hypergeometric_P(b,a+b,c+d,b+d);
132    else if (c<a && c<b && c<d)
133      return  statistics::cdf_hypergeometric_P(c,c+d,a+b,a+c);
134
135    return statistics::cdf_hypergeometric_P(d,c+d,a+b,b+d);
136
137  }
138
139
140  double Fisher::score(const classifier::Target& target, 
141                       const utility::vector& value)
142  {
143    weighted_=false;
144    a_=b_=c_=d_=0;
145    for (size_t i=0; i<target.size(); i++)
146      if (target.binary(i))
147        if (value(i)>value_cutoff_)
148          a_++;
149        else
150          c_++;
151      else
152        if (value(i)>value_cutoff_)
153          b_++;
154        else
155          d_++;
156       
157    return oddsratio(a_,b_,c_,d_); 
158  }
159 
160  double Fisher::score(const classifier::Target& target, 
161                       const classifier::DataLookupWeighted1D& value) 
162  {
163    weighted_=true;
164    a_=b_=c_=d_=0;
165    for (size_t i=0; i<target.size(); i++)
166      if (target.binary(i))
167        if (value.data(i)>value_cutoff_)
168          a_+=value.weight(i);
169        else
170          c_+=value.weight(i);
171      else
172        if (value.data(i)>value_cutoff_)
173          b_+=value.weight(i);
174        else
175          d_+=value.weight(i);
176       
177    return oddsratio(a_,b_,c_,d_); 
178  }
179 
180  double Fisher::score(const classifier::Target& target, 
181                       const utility::vector& value, 
182                       const utility::vector& weight) 
183  {
184    weighted_=true;
185    a_=b_=c_=d_=0;
186    for (size_t i=0; i<target.size(); i++)
187      if (target.binary(i))
188        if (value(i)>value_cutoff_)
189          a_+=weight(i);
190        else
191          c_+=weight(i);
192      else
193        if (value(i)>value_cutoff_)
194          b_+=weight(i);
195        else
196          d_+=weight(i);
197       
198    return oddsratio(a_,b_,c_,d_); 
199  }
200 
201  double Fisher::score(const u_int a, const u_int b, 
202                       const u_int c, const u_int d) 
203  {
204    a_=a;
205    b_=b;
206    c_=c;
207    d_=d;
208    return oddsratio(a,b,c,d); 
209  }
210
211
212  double& Fisher::value_cutoff(void)
213  {
214    return value_cutoff_;
215  }
216
217}}} // of namespace statistics, yat, and theplu
Note: See TracBrowser for help on using the repository browser.