Ignore:
Timestamp:
Jul 12, 2018, 2:43:25 AM (4 years ago)
Author:
Peter
Message:

merge release 0.15.1 into trunk

Location:
trunk
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • trunk

  • trunk/yat/statistics/Fisher.cc

    r3455 r3743  
    55  Copyright (C) 2005 Peter Johansson
    66  Copyright (C) 2006, 2007, 2008 Jari Häkkinen, Peter Johansson
    7   Copyright (C) 2009, 2010, 2011, 2012, 2013, 2014, 2015 Peter Johansson
     7  Copyright (C) 2009, 2010, 2011, 2012, 2013, 2014, 2015, 2018 Peter Johansson
    88
    99  This file is part of the yat library, http://dev.thep.lu.se/yat
     
    8787
    8888  /*
     89   (n)       n!
     90   ( )  = --------
     91   (k)    k!(n-k)!
     92
     93
     94
    8995    ( n )
    90     (k+1)   (k+1)!(n-k)!
    91     ----- = ------------
    92      (n)    k!(n-k-1)!
     96    (k+1)   k!(n-k)!
     97    ----- = ------------ = (n-k) / (k+1)
     98     (n)    (k+1)!(n-k-1)!
    9399     (k)
    94100   */
    95101  double Fisher::choose_ratio(unsigned int n, unsigned int k) const
    96102  {
    97     if (k==0)
    98       return n;
    99     if (k==n-1)
    100       return k+1;
     103    assert(k+1 <= n);
    101104    return static_cast<double>(n-k) / (k+1);
    102105  }
     
    106109                                      unsigned int n2, unsigned int t) const
    107110  {
     111    // P(X = k+1) / P(X = k) =
     112    // (choose(n1, k+1) * choose(n2, t-k-1) / choose(n1+n2, t) ) /
     113    // (choose(n1, k) * choose(n2, t-k) / choose(n1+n2, t) ) =
     114    // choose(n1, k+1) / choose(n1, k) /( choose(n2, t-k)/choose(n2,t-k-1) ) =
     115    // choose(n1,k+1)/choose(n1,k) / (choose(n2,(t-k-1)+1)/choose(n2,t-k-1))
     116    // choose_ratio(n1,k) / choose_ratio(n2, t-k-1)
     117    // where choose_ratio(a, b) = choose(a, b+1) / choose(a,b)
    108118    return choose_ratio(n1, k) / choose_ratio(n2, t-k-1);
    109119  }
     
    240250  {
    241251    // The hypergeometric distribution is unimodal (only one peak) and the
    242     // peak (mode) occurs at: (t+1)*(n1+1)/(n+1)
    243 
    244     double mode = static_cast<double>(a_+c_+1)*static_cast<double>(a_+b_+1)/
    245       static_cast<double>(a_+b_+c_+d_+1);
     252    // peak (mode) occurs at: floor((t+1)*(n1+1)/(n+1))
     253
     254    double mode = std::floor((a_+c_+1.0)*(a_+b_+1.0)/(a_+b_+c_+d_+1.0));
    246255    if (a_ >= mode)
    247256      return p_value_exact(a_, a_+b_, c_+d_, a_+c_);
     
    267276    assert(n1);
    268277    assert(n2);
     278
     279    // special case when k=0 and mode (peak of distribution) is at
     280    // zero as well. If mode is not zero 2x2 table should have been
     281    // mirrored before calling this function.
     282    if (k == 0) {
     283      assert(std::floor((a_+c_+1.0)*(a_+b_+1.0)/(a_+b_+c_+d_+1.0)) == 0.0);
     284      return 1.0;
     285    }
    269286
    270287    assert(k);
Note: See TracChangeset for help on using the changeset viewer.