Changeset 3743 for trunk/yat/statistics/Fisher.cc

Timestamp:

Jul 12, 2018, 2:43:25 AM (4 years ago)

Author:

Peter

Message:

merge release 0.15.1 into trunk

Location:

trunk

Files:

• . (modified) (1 prop)
• yat/statistics/Fisher.cc (modified) (5 diffs)

trunk/yat/statistics/Fisher.cc

@@ -5 +5 @@
   Copyright (C) 2005 Peter Johansson
   Copyright (C) 2006, 2007, 2008 Jari Häkkinen, Peter Johansson
-  Copyright (C) 2009, 2010, 2011, 2012, 2013, 2014, 2015 Peter Johansson
+  Copyright (C) 2009, 2010, 2011, 2012, 2013, 2014, 2015, 2018 Peter Johansson
 
   This file is part of the yat library, http://dev.thep.lu.se/yat
@@ -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);