Changeset 1624 for branches/0.4-stable/yat/statistics/Fisher.cc

Timestamp:

Nov 12, 2008, 11:10:52 PM (15 years ago)

Author:

Peter

Message:

fixes #461. Also modified implementation of cdf_hypergeometric_P, which may cause conflict with modifications done in trunk (refs #87). If so, go with the trunk version (which uses GSL 1.8).

File:

• branches/0.4-stable/yat/statistics/Fisher.cc (modified) (4 diffs)

branches/0.4-stable/yat/statistics/Fisher.cc

@@ -30 +30 @@
 #include <gsl/gsl_cdf.h>
 #include <gsl/gsl_randist.h>
+
+#include <algorithm>
 
 namespace theplu {
@@ -105 +107 @@
   double Fisher::p_value() const
   {
-    if (oddsratio_>1.0)
-      return 2*p_value_one_sided();
-    if (oddsratio_<1.0){
-      // If oddsratio is less than unity
-      // Two-sided p-value is
-      // P(X <= oddsratio) + P(X >= 1/oddsratio) = 
-      // P(X <= oddsratio) + P(X <= oddsratio)
-      // 2 * (P(X < oddsratio) + P(X = oddsratio))
-      // 2 * (1 - P(X >= oddsratio) + P(X = oddsratio))
-      if (calculate_p_exact())
-        return 2*(1-p_value_one_sided()+
-                  gsl_ran_hypergeometric_pdf(a_,a_+b_,c_+d_, a_+c_)
-                  );
-      // if p is calculated approximatively correction is not needed
-      return 2*(1-p_value_one_sided());
-
-    }
-    return 1.0;
+    if (calculate_p_exact())
+      return p_value_exact();
+    return p_value_approximative();
   }
 
@@ -128 +115 @@
   double Fisher::p_value_one_sided() const
   {
-    if ( calculate_p_exact() )
-      return p_value_exact();
-    return p_value_approximative();
+    if (!calculate_p_exact()) {
+      if (oddsratio_<1.0)
+        return 1.0-p_value_approximative()/2.0;
+      return p_value_approximative()/2.0;
+    }
+    return statistics::cdf_hypergeometric_P(c_, c_+d_, a_+b_, a_+c_);
   }
 
@@ -141 +131 @@
   double Fisher::p_value_exact() const
   { 
-    // Since the calculation is symmetric and cdf_hypergeometric_P
-    // loops up to k we choose the smallest number to be k and mirror
-    // the matrix. This choice makes the p-value two-sided.
+    double res=0;
+    double a, c, tmp;
+    expected(a, tmp, c, tmp);
+    // sum P(x) over x for which abs(x-E(a))>=abs(a-E(a))  
 
-    if (a_<b_ && a_<c_ && a_<d_)
-      return statistics::cdf_hypergeometric_P(a_,a_+b_,c_+d_,a_+c_);
-    else if (b_<a_ && b_<c_ && b_<d_)
-      return  statistics::cdf_hypergeometric_P(b_,a_+b_,c_+d_,b_+d_);
-    else if (c_<a_ && c_<b_ && c_<d_)
-      return  statistics::cdf_hypergeometric_P(c_,c_+d_,a_+b_,a_+c_);
+    // counting left tail
+    int k = static_cast<int>(a - std::abs(a_-a));
+    if (k>=0)
+      res = cdf_hypergeometric_P(k, a_+b_, c_+d_, a_+c_);
+    // counting right tail
+    k = static_cast<int>(c - std::abs(c_-c));
+    if (k>=0)
+      res += cdf_hypergeometric_P(k, c_+d_, a_+b_, a_+c_);
 
-    return statistics::cdf_hypergeometric_P(d_,c_+d_,a_+b_,b_+d_);
-
+    // avoid p>1 due to double counting middle one
+    return std::min(1.0, res);
   }