Changeset 1624 for branches

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).

Location:

branches/0.4-stable

Files:

• build_support/version.m4 (modified) (1 diff)
• test/fisher_test.cc (modified) (2 diffs)
• yat/statistics/Fisher.cc (modified) (4 diffs)
• yat/statistics/Fisher.h (modified) (1 diff)
• yat/statistics/utility.cc (modified) (2 diffs)

branches/0.4-stable/build_support/version.m4

@@ -54 +54 @@
 # yat-0.4.1  1:0:0
 # yat-0.4.2  1:1:0
+# yat-0.4.3  1:2:0
 #
-m4_define([LT_VERSION], [1:1:0])
+m4_define([LT_VERSION], [1:2:0])
 
 

branches/0.4-stable/test/fisher_test.cc

@@ -24 +24 @@
 #include "yat/statistics/Fisher.h"
 
+using namespace theplu::yat;
+void test_p_value(test::Suite&);
+void test_p_value_approximative(test::Suite&);
+void test_p_value_exact(test::Suite&);
+
 int main(int argc, char* argv[])
 {  
-  using namespace theplu::yat;
   test::Suite suite(argc, argv);
 
@@ -60 +64 @@
     suite.err() << "minimum_size failed\n";
   }
+  test_p_value(suite);
   return suite.return_value();
 }
+
+
+void test_p_value(test::Suite& suite)
+{
+  test_p_value_exact(suite);
+  test_p_value_approximative(suite);
+}
+
+
+void test_p_value_approximative(test::Suite& suite)
+{
+  suite.err() << "testing p_value_approximative\n";
+  statistics::Fisher f;
+  f.minimum_size() = 0;
+  f.oddsratio(10,20,20,50);
+  suite.add(suite.equal(f.p_value(), 2*f.p_value_one_sided()));
+  f.oddsratio(10,20,10,20);
+  suite.add(suite.equal(f.p_value(), 1.0));
+  suite.add(suite.equal(f.p_value_one_sided(), 0.5));
+}
+
+
+void test_p_value_exact(test::Suite& suite)
+{
+  suite.err() << "test p_value_exact\n";
+  statistics::Fisher f;
+  f.minimum_size() = 1000;
+  f.oddsratio(10,20,10,20);
+  suite.add(suite.equal(f.p_value(), 1.0));
+
+  f.oddsratio(10, 20, 20, 200);
+  suite.add(suite.equal(f.p_value(), 0.0008119060627676223));
+  suite.add(suite.equal(f.p_value_one_sided(), 0.0008119060627676223));
+
+  // testing symmetry
+  statistics::Fisher f2;
+  f2.minimum_size() = 1000;
+  f2.oddsratio(20, 200, 10, 20);
+  suite.add(suite.equal(f2.p_value(), f.p_value()));
+
+  f.oddsratio(1, 1, 1, 2);
+  suite.add(suite.equal(f.p_value(), 1.0));
+  suite.add(suite.equal(f.p_value_one_sided(), 0.7));
+
+  f.oddsratio(1, 1, 2, 1);
+  suite.add(suite.equal(f.p_value(), 1.0));
+  suite.add(suite.equal(f.p_value_one_sided(), 0.9));
+
+}

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);
   }
 

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

@@ -118 +118 @@
     const unsigned int& minimum_size(void) const;
 
-    ///
-    /// If oddsratio is larger than unity, two-sided p-value is equal
-    /// to 2*p_value_one_sided(). If oddsratio is smaller than unity
-    /// two-sided p-value is equal to 2*(1-p_value_one_sided()). If
-    /// oddsratio is unity two-sided p-value is equal to unity.
-    ///
-    /// If all elements in table is at least minimum_size(), a Chi2
-    /// approximation is used.
-    ///
-    /// @return 2-sided p-value
-    ///
+    /**
+       If all elements in table is at least minimum_size(), a Chi2
+       approximation is used.
+       
+       Otherwise a two-sided p-value is calculated using the
+       hypergeometric distribution
+       \f$ \sum_k P(k) \f$ where summation runs over \a k such that
+       \f$ |k-<a>| \ge |a-<a>| \f$.
+
+       \return two-sided p-value
+    */
     double p_value() const;
     

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

@@ -31 +31 @@
 #include <gsl/gsl_statistics_double.h>
 
+#include <algorithm>
 #include <cassert>
 
@@ -41 +42 @@
   {
     double p=0;
-    for (size_t i=0; i<=k; i++)
+    size_t top = std::min(k, std::min(n1, t));
+    for (size_t i=std::max(0, static_cast<int>(t-n2)); i<=top; i++)
       p+= gsl_ran_hypergeometric_pdf(i, n1, n2, t);
     return p;