Changeset 449
 Timestamp:
 Dec 15, 2005, 9:06:10 PM (17 years ago)
 Location:
 trunk/lib/statistics
 Files:

 2 edited
Legend:
 Unmodified
 Added
 Removed

trunk/lib/statistics/Fisher.cc
r448 r449 5 5 #include <c++_tools/statistics/utility.h> 6 6 7 #include <gsl/gsl_cdf.h> 8 #include <gsl/gsl_randist.h> 9 7 10 namespace theplu { 8 11 namespace statistics { … … 10 13 11 14 Fisher::Fisher(bool b) 12 : Score(b), a_(0), b_(0), c_(0), d_(0), cutoff_column_(0), cutoff_row_(0) 15 : Score(b), a_(0), b_(0), c_(0), d_(0), cutoff_column_(0), cutoff_row_(0), 16 oddsratio_(1.0) 13 17 { 14 18 } 15 19 16 20 17 void Fisher::expected(double& a, double& b, double& c, double& d) 21 double Fisher::Chi2() const 22 { 23 double a,b,c,d; 24 expected(a,b,c,d); 25 return (aa_)*(aa_)/a + (bb_)*(bb_)/b + 26 (cc_)*(cc_)/c + (dd_)*(dd_)/d; 27 } 28 29 30 void Fisher::expected(double& a, double& b, double& c, double& d) const 18 31 { 19 32 double N = a_+b_+c_+d_; … … 27 40 const double b, 28 41 const double c, 29 const double d) const42 const double d) 30 43 { 31 44 // If a column sum or a row sum is zero, the table is nonsense … … 33 46 //Peter, should through exception 34 47 std::cerr << "Warning: Fisher: Table is not valid\n"; 35 return 1.0;48 return oddsratio_ = 1.0; 36 49 } 37 50 38 double ratio=(a*d)/(b*d);39 if (absolute_ && ratio<1)40 return 1/ ratio;41 return ratio;51 oddsratio_=(a*d)/(b*d); 52 if (absolute_ && oddsratio_<1) 53 return 1/oddsratio_; 54 return oddsratio_; 42 55 } 43 56 … … 54 67 p=p_value_approximative(); 55 68 56 if (absolute_) 57 return fabs(2*(p.5)); 69 if (!absolute_){ 70 p=p/2; 71 if (oddsratio_<0.5){ 72 // last term because >= not equal to !(<=) 73 return 1p+gsl_ran_hypergeometric_pdf(a_, a_+b_, c_+d_, a_+c_); 74 } 75 } 58 76 return p; 59 77 } … … 62 80 double Fisher::p_value_approximative() const 63 81 { 64 return 1.0;82 return gsl_cdf_chisq_Q(Chi2(), 1.0); 65 83 } 66 67 84 68 85 double Fisher::p_value_exact() const … … 70 87 // Since the calculation is symmetric and cdf_hypergeometric_P 71 88 // loops up to k we choose the smallest number to be k and mirror 72 // the matrix. 89 // the matrix. This choice makes the pvalue twosided. 73 90 if (a_<b_ && a_<c_ && a_<d_) 74 91 return statistics::cdf_hypergeometric_P(a_,a_+b_,c_+d_,a_+c_); 75 92 else if (b_<a_ && b_<c_ && b_<d_) 76 return statistics::cdf_hypergeometric_P(b_,a_+b_,c_+d_,b_+d_);93 return statistics::cdf_hypergeometric_P(b_,a_+b_,c_+d_,b_+d_); 77 94 else if (c_<a_ && c_<b_ && c_<d_) 78 return statistics::cdf_hypergeometric_P(c_,c_+d_,a_+b_,a_+c_); 79 else 80 return statistics::cdf_hypergeometric_P(d_,c_+d_,a_+b_,b_+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 81 99 } 82 100 
trunk/lib/statistics/Fisher.h
r448 r449 58 58 59 59 /// 60 /// @return Chi2 score 61 /// 62 double Chi2(void) const; 63 64 /// 60 65 /// Cutoff sets the limit whether a value should go into the left 61 66 /// or the right column. @see score … … 77 82 /// a' = (a+c)(a+b)/(a+b+c+d) 78 83 /// 79 void expected(double& a, double& b, double& c, double& d) ;84 void expected(double& a, double& b, double& c, double& d) const; 80 85 81 86 /// … … 136 141 private: 137 142 double oddsratio(const double a, const double b, 138 const double c, const double d) const;143 const double c, const double d); 139 144 145 // twosided 140 146 double p_value_approximative(void) const; 147 //twosided 141 148 double p_value_exact(void) const; 142 149 … … 150 157 double cutoff_row_; 151 158 u_int minimum_size_; 152 159 double oddsratio_; 153 160 }; 154 161
Note: See TracChangeset
for help on using the changeset viewer.