source: trunk/c++_tools/statistics/Fisher.cc @ 623

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

fixes #112 and refs #123 added overloaded function score taking Target and DataLookupWeighted1D, which is needed for InputRanker?.

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 4.7 KB
Line 
1// $Id: Fisher.cc 623 2006-09-05 02:13:12Z peter $
2
3#include <c++_tools/statistics/Fisher.h>
4#include <c++_tools/statistics/Score.h>
5#include <c++_tools/statistics/utility.h>
6#include <c++_tools/classifier/DataLookupWeighted1D.h>
7#include <c++_tools/classifier/Target.h>
8
9#include <gsl/gsl_cdf.h>
10#include <gsl/gsl_randist.h>
11
12namespace theplu {
13namespace statistics {
14
15 
16  Fisher::Fisher(bool b)
17    : Score(b), a_(0), b_(0), c_(0), d_(0), oddsratio_(1.0), value_cutoff_(0)
18  {
19  }
20
21
22  double Fisher::Chi2() const
23  {
24    double a,b,c,d;
25    expected(a,b,c,d);
26    return (a-a_)*(a-a_)/a + (b-b_)*(b-b_)/b + 
27      (c-c_)*(c-c_)/c + (d-d_)*(d-d_)/d;
28  }
29
30
31  void Fisher::expected(double& a, double& b, double& c, double& d) const
32  {
33    double N = a_+b_+c_+d_;
34    a =((a_+b_)*(a_+c_)) / N;
35    b =((a_+b_)*(b_+d_)) / N;
36    c =((c_+d_)*(a_+c_)) / N;
37    d =((c_+d_)*(b_+d_)) / N;
38  }
39
40  double Fisher::oddsratio(const double a,
41                           const double b,
42                           const double c,
43                           const double d) 
44  {
45    // If a column sum or a row sum is zero, the table is nonsense
46    if ((a==0 || d==0) && (c==0 || b==0)){
47      //Peter, should throw exception
48      std::cerr << "Warning: Fisher: Table is not valid\n";
49      return oddsratio_ = 1.0;
50    }
51
52    oddsratio_=(a*d)/(b*d);
53    if (absolute_ && oddsratio_<1)
54      return 1/oddsratio_;
55    return oddsratio_;
56  }
57
58
59  double Fisher::p_value() const
60  {
61    double p=1;
62    if ( ( a_<minimum_size_ || b_<minimum_size_ || 
63           c_<minimum_size_ || d_<minimum_size_) && !weighted_ ){
64
65      p=p_value_exact();
66    }
67    else
68      p=p_value_approximative();
69
70    if (!absolute_){
71      p=p/2;
72      if (oddsratio_<0.5){
73        // last term because >= not equal to !(<=)
74        u_int a = static_cast<u_int>(a_);
75        u_int b = static_cast<u_int>(b_);
76        u_int c = static_cast<u_int>(c_);
77        u_int d = static_cast<u_int>(d_);
78        return 1-p+gsl_ran_hypergeometric_pdf(a, a+b, c+d, a+c);
79      }
80    }
81    return p;
82  }
83
84
85  double Fisher::p_value_approximative() const
86  {
87    return gsl_cdf_chisq_Q(Chi2(), 1.0);
88  }
89
90  double Fisher::p_value_exact() const
91  { 
92    u_int a = static_cast<u_int>(a_);
93    u_int b = static_cast<u_int>(b_);
94    u_int c = static_cast<u_int>(c_);
95    u_int d = static_cast<u_int>(d_);
96   
97    // Since the calculation is symmetric and cdf_hypergeometric_P
98    // loops up to k we choose the smallest number to be k and mirror
99    // the matrix. This choice makes the p-value two-sided.
100
101    if (a<b && a<c && a<d)
102      return statistics::cdf_hypergeometric_P(a,a+b,c+d,a+c);
103    else if (b<a && b<c && b<d)
104      return  statistics::cdf_hypergeometric_P(b,a+b,c+d,b+d);
105    else if (c<a && c<b && c<d)
106      return  statistics::cdf_hypergeometric_P(c,c+d,a+b,a+c);
107
108    return statistics::cdf_hypergeometric_P(d,c+d,a+b,b+d);
109
110  }
111
112
113    double Fisher::score(const classifier::Target& target, 
114                         const utility::vector& value)
115  {
116    weighted_=false;
117    a_=b_=c_=d_=0;
118    for (size_t i=0; i<target.size(); i++)
119      if (target.binary(i))
120        if (value(i)>value_cutoff_)
121          a_++;
122        else
123          c_++;
124      else
125        if (value(i)>value_cutoff_)
126          b_++;
127        else
128          d_++;
129       
130    // If a column sum or a row sum is zero, the table is non-sense
131    if ((a_==0 || d_==0) && (c_==0 || b_==0)){
132      std::cerr << "Warning: Fisher: Table is not valid\n";
133      return 1;
134    }
135
136    return oddsratio(a_,b_,c_,d_); 
137  }
138 
139  double Fisher::score(const classifier::Target& target, 
140                       const classifier::DataLookupWeighted1D& value) 
141  {
142    weighted_=true;
143    a_=b_=c_=d_=0;
144    for (size_t i=0; i<target.size(); i++)
145      if (target.binary(i))
146        if (value.data(i)>value_cutoff_)
147          a_+=value.weight(i);
148        else
149          c_+=value.weight(i);
150      else
151        if (value.data(i)>value_cutoff_)
152          b_+=value.weight(i);
153        else
154          d_+=value.weight(i);
155       
156    // If a column sum or a row sum is zero, the table is non-sense
157    if ((a_==0 || d_==0) && (c_==0 || b_==0)){
158      // Peter should throw an exception here
159      std::cerr << "Warning: Fisher: Table is not valid\n";
160      return 1;
161    }
162
163    return oddsratio(a_,b_,c_,d_); 
164  }
165 
166    double Fisher::score(const classifier::Target& target, 
167                         const utility::vector& value, 
168                         const utility::vector& weight) 
169  {
170    weighted_=true;
171    a_=b_=c_=d_=0;
172    for (size_t i=0; i<target.size(); i++)
173      if (target.binary(i))
174        if (value(i)>value_cutoff_)
175          a_+=weight(i);
176        else
177          c_+=weight(i);
178      else
179        if (value(i)>value_cutoff_)
180          b_+=weight(i);
181        else
182          d_+=weight(i);
183       
184    // If a column sum or a row sum is zero, the table is non-sense
185    if ((a_==0 || d_==0) && (c_==0 || b_==0)){
186      // Peter should throw an exception
187      std::cerr << "Warning: Fisher: Table is not valid\n";
188      return 1;
189    }
190
191    return oddsratio(a_,b_,c_,d_); 
192  }
193 
194  double Fisher::score(const u_int a, const u_int b, 
195                       const u_int c, const u_int d) 
196  {
197    a_=a;
198    b_=b;
199    c_=c;
200    d_=d;
201    return oddsratio(a,b,c,d); 
202  }
203
204}} // of namespace statistics and namespace theplu
Note: See TracBrowser for help on using the repository browser.