source: trunk/lib/statistics/Fisher.cc @ 449

Last change on this file since 449 was 449, checked in by Peter, 16 years ago

added approximative p-value calculation using Chi2, fixed some bugs, and extended the interface slightly

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 4.0 KB
Line 
1// $Id: Fisher.cc 449 2005-12-15 20:06:10Z peter $
2
3#include <c++_tools/statistics/Fisher.h>
4#include <c++_tools/statistics/Score.h>
5#include <c++_tools/statistics/utility.h>
6
7#include <gsl/gsl_cdf.h>
8#include <gsl/gsl_randist.h>
9
10namespace theplu {
11namespace statistics {
12
13 
14  Fisher::Fisher(bool b)
15    : Score(b), a_(0), b_(0), c_(0), d_(0), cutoff_column_(0), cutoff_row_(0), 
16      oddsratio_(1.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 through 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_){
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        return 1-p+gsl_ran_hypergeometric_pdf(a_, a_+b_, c_+d_, a_+c_);
74      }
75    }
76    return p;
77  }
78
79
80  double Fisher::p_value_approximative() const
81  {
82    return gsl_cdf_chisq_Q(Chi2(), 1.0);
83  }
84
85  double Fisher::p_value_exact() const
86  {
87    // Since the calculation is symmetric and cdf_hypergeometric_P
88    // loops up to k we choose the smallest number to be k and mirror
89    // the matrix. This choice makes the p-value two-sided.
90    if (a_<b_ && a_<c_ && a_<d_)
91      return statistics::cdf_hypergeometric_P(a_,a_+b_,c_+d_,a_+c_);
92    else if (b_<a_ && b_<c_ && b_<d_)
93      return  statistics::cdf_hypergeometric_P(b_,a_+b_,c_+d_,b_+d_);
94    else if (c_<a_ && c_<b_ && c_<d_)
95      return  statistics::cdf_hypergeometric_P(c_,c_+d_,a_+b_,a_+c_);
96
97    return statistics::cdf_hypergeometric_P(d_,c_+d_,a_+b_,b_+d_);
98
99  }
100
101
102  double Fisher::score(const gslapi::vector& x, 
103                       const gslapi::vector& y, 
104                       const std::vector<size_t>& train_set)
105  {
106    if (!train_set.size()){
107      train_set_.resize(0);
108      for (size_t i=0; i<target_.size(); i++)
109        train_set_.push_back(i); 
110    }
111    else
112      train_set_ = train_set;
113    a_=b_=c_=d_=0;
114    for (size_t i=0; i<train_set_.size(); i++)
115      if (x(train_set_[i])<cutoff_column_)
116        if (y(train_set_[i])<cutoff_row_)
117          a_++;
118        else
119          c_++;
120      else
121        if (y(train_set_[i])<cutoff_row_)
122          b_++;
123        else
124          d_++;
125       
126    // If a column sum or a row sum is zero, the table is nonsense
127    if ((a_==0 || d_==0) && (c_==0 || b_==0)){
128      std::cerr << "Warning: Fisher: Table is not valid\n";
129      return 0;
130    }
131
132    return oddsratio(a_,b_,c_,d_); 
133  }
134 
135  double Fisher::score(const gslapi::vector& x, 
136                       const gslapi::vector& y,             
137                       const gslapi::vector& w, 
138                       const std::vector<size_t>& train_set)
139  {
140    if (!train_set.size()){
141      train_set_.resize(0);
142      for (size_t i=0; i<target_.size(); i++)
143        train_set_.push_back(i); 
144    }
145    else
146      train_set_ = train_set;
147
148    double a=0;
149    double b=0;
150    double c=0;
151    double d=0;
152
153    for (size_t i=0; i<train_set_.size(); i++)
154      if (x(train_set_[i]) < cutoff_column_)
155        if (y(train_set_[i]) < cutoff_row_)
156          a+=w(train_set_[i]);
157        else
158          c+=w(train_set_[i]);
159      else
160        if (y(train_set_[i]) < cutoff_column_)
161          b+=w(train_set_[i]);
162        else
163          d+=w(train_set_[i]);
164
165    return oddsratio(a_,b_,c_,d_); 
166  }
167
168  double Fisher::score(const u_int a, const u_int b, 
169                       const u_int c, const u_int d) 
170  {
171    a_=a;
172    b_=b;
173    c_=c;
174    d_=d;
175    return oddsratio(a,b,c,d); 
176  }
177
178}} // of namespace statistics and namespace theplu
Note: See TracBrowser for help on using the repository browser.