Changeset 3739 for branches/0.15-stable

Timestamp:

Jul 5, 2018, 7:59:40 AM (5 years ago)

Author:

Peter

Message:

fixes #908

Location:

branches/0.15-stable

Files:

• test/fisher.cc (modified) (4 diffs)
• yat/statistics/Fisher.cc (modified) (4 diffs)

branches/0.15-stable/test/fisher.cc

@@ -26 +26 @@
 #include "yat/statistics/Fisher.h"
 
+#include <gsl/gsl_cdf.h>
+#include <gsl/gsl_randist.h>
 #include <climits>
 
@@ -34 +36 @@
 void test_large_numbers(test::Suite&);
 void test_yates(test::Suite&);
+void test_cdf_hypergeometric_Q(test::Suite&);
+void test_ticket908(test::Suite&);
 
 int main(int argc, char* argv[])
@@ -73 +77 @@
   test_large_numbers(suite);
   test_yates(suite);
+  test_cdf_hypergeometric_Q(suite);
+  test_ticket908(suite);
   return suite.return_value();
 }
@@ -168 +174 @@
   }
 }
+
+
+void test_cdf_hypergeometric_Q(test::Suite& suite)
+{
+  suite.out() << "\ntest " << __func__ << "\n";
+  unsigned int n1 = 16;
+  unsigned int n2 = 3;
+  unsigned int t = 16;
+  for (unsigned int k=0; k<=t; ++k) {
+    // P(X=k)
+    double p = gsl_ran_hypergeometric_pdf(k, n1, n2, t);
+    // P(X<=k)
+    double P = gsl_cdf_hypergeometric_P(k, n1, n2, t);
+    // P(X>k)
+    double Q = gsl_cdf_hypergeometric_Q(k, n1, n2, t);
+    suite.out() << k << " " << p << " " << P << " " << Q << "\n";
+    double sum = 0;
+    for (unsigned int x=0; x<=k; ++x)
+      sum += gsl_ran_hypergeometric_pdf(x, n1, n2, t);
+    if (!suite.equal(sum, P, 20)) {
+      suite.add(false);
+      suite.err() << "error: gsl_cdf_hypergeometric_P for k=" << k << "\n";
+    }
+    if (!suite.equal(Q, 1-P, 1000)) {
+      suite.add(false);
+      suite.err() << "failed: Q = 1-P for k=" << k << "\n";
+    }
+  }
+}
+
+
+// test bug #908 reported in http://dev.thep.lu.se/yat/ticket/908
+void test_ticket908(test::Suite& suite)
+{
+  suite.out() << "===\ntesting ticket 908\n";
+  statistics::Fisher f;
+  int a1 = 13;
+  int a2 = 3;
+  int b1 = 3;
+  int b2 = 0;
+
+  double correct_p = 0;
+  unsigned int n1 = a1 + a2;
+  unsigned int n2 = b1 + b2;
+  unsigned int t = a1+b1;
+  for (unsigned int k=0; k<=n1; ++k) {
+    double P = gsl_ran_hypergeometric_pdf(k, n1, n2, t);
+    if (P <=  gsl_ran_hypergeometric_pdf(a1, n1, n2, t)) {
+      suite.out() << "counting " << k << " " << P << "\n";
+      correct_p += P;
+    }
+  }
+  statistics::Fisher fisher;
+  std::vector<double> p;
+  fisher.oddsratio(a1, a2, b1, b2);
+  p.push_back(fisher.p_value());
+  fisher.oddsratio(a2, a1, b2, b1);
+  p.push_back(fisher.p_value());
+  fisher.oddsratio(b1, b2, a1, a2);
+  p.push_back(fisher.p_value());
+  fisher.oddsratio(b2, b1, a2, a1);
+  p.push_back(fisher.p_value());
+  for (size_t i=0; i<p.size(); ++i) {
+    if (!suite.add(suite.equal(correct_p, p[i], 20)))
+      suite.err() << "error: incorrect p: " << p[i] << " for case "
+                  << i << "; expected: " << correct_p << "\n";
+  }
+}

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

@@ -87 +87 @@
 
   /*
+   (n)       n!
+   ( )  = --------
+   (k)    k!(n-k)!
+
+
+
     ( n )
-    (k+1)   (k+1)!(n-k)!
-    ----- = ------------
-     (n)    k!(n-k-1)!
+    (k+1)   k!(n-k)!
+    ----- = ------------ = (n-k) / (k+1)
+     (n)    (k+1)!(n-k-1)!
      (k)
    */
   double Fisher::choose_ratio(unsigned int n, unsigned int k) const
   {
-    if (k==0)
-      return n;
-    if (k==n-1)
-      return k+1;
+    assert(k+1 <= n);
     return static_cast<double>(n-k) / (k+1);
   }
@@ -106 +109 @@
                                       unsigned int n2, unsigned int t) const
   {
+    // P(X = k+1) / P(X = k) =
+    // (choose(n1, k+1) * choose(n2, t-k-1) / choose(n1+n2, t) ) /
+    // (choose(n1, k) * choose(n2, t-k) / choose(n1+n2, t) ) =
+    // choose(n1, k+1) / choose(n1, k) /( choose(n2, t-k)/choose(n2,t-k-1) ) =
+    // choose(n1,k+1)/choose(n1,k) / (choose(n2,(t-k-1)+1)/choose(n2,t-k-1))
+    // choose_ratio(n1,k) / choose_ratio(n2, t-k-1)
+    // where choose_ratio(a, b) = choose(a, b+1) / choose(a,b)
     return choose_ratio(n1, k) / choose_ratio(n2, t-k-1);
   }
@@ -240 +250 @@
   {
     // The hypergeometric distribution is unimodal (only one peak) and the
-    // peak (mode) occurs at: (t+1)*(n1+1)/(n+1)
-
-    double mode = static_cast<double>(a_+c_+1)*static_cast<double>(a_+b_+1)/
-      static_cast<double>(a_+b_+c_+d_+1);
+    // peak (mode) occurs at: floor((t+1)*(n1+1)/(n+1))
+
+    double mode = std::floor((a_+c_+1.0)*(a_+b_+1.0)/(a_+b_+c_+d_+1.0));
     if (a_ >= mode)
       return p_value_exact(a_, a_+b_, c_+d_, a_+c_);
@@ -267 +276 @@
     assert(n1);
     assert(n2);
+
+    // special case when k=0 and mode (peak of distribution) is at
+    // zero as well. If mode is not zero 2x2 table should have been
+    // mirrored before calling this function.
+    if (k == 0) {
+      assert(std::floor((a_+c_+1.0)*(a_+b_+1.0)/(a_+b_+c_+d_+1.0)) == 0.0);
+      return 1.0;
+    }
 
     assert(k);