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

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

Removed gslapi namespace and put the code into utility namespace.
Removed unneccesary #includes, and added needed #includes.

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