Changeset 3441

Timestamp:

Nov 23, 2015, 9:20:08 AM (7 years ago)

Author:

Peter

Message:

closes #841

Location:

trunk

Files:

• test/fisher.cc (modified) (1 diff)
• yat/statistics/Fisher.cc (modified) (4 diffs)
• yat/statistics/Fisher.h (modified) (2 diffs)

trunk/test/fisher.cc

@@ -126 +126 @@
   f.minimum_size() = 1000;
   f.oddsratio(10,20,10,20);
-  suite.add(suite.equal(f.p_value(), 1.0));
+  suite.add(suite.equal(f.p_value(), 1.0, 20));
 
   f.oddsratio(10, 20, 20, 200);

trunk/yat/statistics/Fisher.cc

@@ -86 +86 @@
 
 
+  /*
+    ( n )
+    (k+1)   (k+1)!(n-k)!
+    ----- = ------------
+     (n)    k!(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;
+    return static_cast<double>(n-k) / (k+1);
+  }
+
+
+  double Fisher::hypergeometric_ratio(unsigned int k, unsigned int n1,
+                                      unsigned int n2, unsigned int t) const
+  {
+    return choose_ratio(n1, k) / choose_ratio(n2, t-k-1);
+  }
+
+  double Fisher::lower_tail(unsigned int k, unsigned int n1, unsigned int n2,
+                            unsigned int t) const
+  {
+    double sum = 0;
+
+    // P(X=k)
+    double p0 = gsl_ran_hypergeometric_pdf(k, n1, n2, t);
+
+    // minimum possible outcome for X is max(0, t-n2)
+    unsigned int i = std::max(n2, t)-n2;
+    double p = gsl_ran_hypergeometric_pdf(i, n1, n2, t);
+    // avoid double dipping P(k)
+    for ( ; i<k; ++i) {
+      if (p<=p0)
+        sum += p;
+      else
+        break;
+
+      // calculate p(i+1) using recursive function
+      p *= hypergeometric_ratio(i, n1, n2, t);
+    }
+
+    return sum;
+  }
+
+
   unsigned int& Fisher::minimum_size(void)
   {
@@ -131 +180 @@
       return p_value_approximative()/2.0;
     }
+    return p_left_exact();
+  }
+
+
+  double Fisher::p_left_exact(void) const
+  {
     // check for overflow
     assert(c_ <= c_+d_ && d_ <= c_+d_);
@@ -155 +210 @@
       return p_value_approximative()/2.0;
     }
+    return p_right_exact();
+  }
+
+
+  double Fisher::p_right_exact(void) const
+  {
     // check for overflow
     assert(c_ <= c_+d_ && d_ <= c_+d_);
@@ -178 +239 @@
   double Fisher::p_value_exact(void) const
   {
-    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))
-
-    // 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_);
-
-    // avoid p>1 due to double counting middle one
-    return std::min(1.0, res);
+    // 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);
+    if (a_ >= mode)
+      return p_value_exact(a_, a_+b_, c_+d_, a_+c_);
+    return p_value_exact(c_, c_+d_, a_+b_, a_+c_);
+  }
+
+
+  /*
+    a | b | n1
+    c | d | n2
+    -------------------
+    t |   | n
+
+
+    The p is calculated as the sum of all P(x) for all x such that
+    P(x)<=P(k)
+  */
+  double Fisher::p_value_exact(unsigned int k, unsigned int n1,
+                               unsigned int n2, unsigned int t) const
+  {
+    assert(k<=t);
+    assert(t<=n1+n2);
+    assert(n1);
+    assert(n2);
+
+    assert(k);
+    // P(X >= k) = P(X > k-1)
+    double sum = gsl_cdf_hypergeometric_Q(k-1, n1, n2, t);
+    // calculate the other tail
+    sum += lower_tail(k, n1, n2, t);
+    return sum;
   }
 

trunk/yat/statistics/Fisher.h

@@ -148 +148 @@
        hypergeometric distribution
        \f$ \sum_k P(k) \f$ where summation runs over \a k such that
-       \f$ |k-<a>| \ge |a-<a>| \f$.
+       \f$ P(k) \le P(a) \f$.
 
        \return two-sided p-value
@@ -195 +195 @@
     // two-sided
     double p_value_approximative(void) const;
-    //two-sided
+    // two-sided
     double p_value_exact(void) const;
+    // calculate two-sided p-value to get k (or fewer) wins when
+    // drawing t samples out of of a population of n1 wins and n2 losses
+    double p_value_exact(unsigned int k, unsigned int n1, unsigned int n2,
+                         unsigned int t) const;
+
+    double lower_tail(unsigned int k, unsigned int n1, unsigned int n2,
+                      unsigned int t) const;
+
+    // return P(X=k+1) / P(X=k)
+    double hypergeometric_ratio(unsigned int k, unsigned int n1,
+                                unsigned int n2, unsigned int t) const;
+    double choose_ratio(unsigned int n, unsigned int k) const;
+
+    double p_left_exact(void) const;
+    double p_right_exact(void) const;
 
     double yates(unsigned int o, double e) const;