source: trunk/yat/statistics/Fisher.cc @ 3743

// $Id: Fisher.cc 3743 2018-07-12 00:43:25Z peter $

/*
  Copyright (C) 2004 Jari Häkkinen, Peter Johansson
  Copyright (C) 2005 Peter Johansson
  Copyright (C) 2006, 2007, 2008 Jari Häkkinen, 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

  The yat library is free software; you can redistribute it and/or
  modify it under the terms of the GNU General Public License as
  published by the Free Software Foundation; either version 3 of the
  License, or (at your option) any later version.

  The yat library is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
  General Public License for more details.

  You should have received a copy of the GNU General Public License
  along with yat. If not, see <http://www.gnu.org/licenses/>.
*/

#include <config.h>

#include "Fisher.h"
#include "utility.h"

#include "yat/utility/Exception.h"

#include <gsl/gsl_cdf.h>
#include <gsl/gsl_randist.h>

#include <algorithm>
#include <cassert>

namespace theplu {
namespace yat {
namespace statistics {


  Fisher::Fisher(bool yates)
    : a_(0), b_(0), c_(0), d_(0), minimum_size_(10), oddsratio_(1.0),
      yates_(yates)
  {
  }


  Fisher::~Fisher()
  {
  }


  bool Fisher::calculate_p_exact(void) const
  {
    return ( a_<minimum_size_ || b_<minimum_size_ ||
             c_<minimum_size_ || d_<minimum_size_);
  }

  double Fisher::Chi2(void) const
  {
    double a,b,c,d;
    expected(a,b,c,d);
    if (yates_)
      return yates(a_, a) + yates(b_, b) + yates(c_, c) + yates(d_, d);

    return (a-a_)*(a-a_)/a + (b-b_)*(b-b_)/b +
      (c-c_)*(c-c_)/c + (d-d_)*(d-d_)/d;
  }


  void Fisher::expected(double& a, double& b, double& c, double& d) const
  {
    // use floting point arithmetic to avoid overflow
    double a1 = a_;
    double b1 = b_;
    double c1 = c_;
    double d1 = d_;
    double N = a1 + b1 + c1 + d1;
    a =((a1+b1)*(a1+c1)) / N;
    b =((a1+b1)*(b1+d1)) / N;
    c =((c1+d1)*(a1+c1)) / N;
    d =((c1+d1)*(b1+d1)) / N;
  }


  /*
   (n)       n!
   ( )  = --------
   (k)    k!(n-k)!



    ( n )
    (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
  {
    assert(k+1 <= n);
    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
  {
    // 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);
  }

  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)
  {
    return minimum_size_;
  }


  const unsigned int& Fisher::minimum_size(void) const
  {
    return minimum_size_;
  }


  double Fisher::oddsratio(const unsigned int a,
                           const unsigned int b,
                           const unsigned int c,
                           const unsigned int d)
  {
    // If a column sum or a row sum is zero, the table is nonsense
    if ((a==0 || d==0) && (c==0 || b==0)){
      throw utility::runtime_error("Table in Fisher is not valid\n");
    }
    a_ = a;
    b_ = b;
    c_ = c;
    d_ = d;

    oddsratio_= (static_cast<double>(a) / static_cast<double>(b) *
                 static_cast<double>(d) / static_cast<double>(c) );
    return oddsratio_;
  }


  double Fisher::oddsratio(void) const
  {
    return oddsratio_;
  }


  double Fisher::p_left(void) const
  {
    if (!calculate_p_exact()) {
      if (oddsratio_>1.0)
        return 1.0-p_value_approximative()/2.0;
      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_);
    assert(a_ <= a_+b_ && b_ <= a_+b_);
    assert(a_ <= a_+c_ && c_ <= a_+c_);

    return cdf_hypergeometric_P(a_, a_+b_, c_+d_, a_+c_);
  }


  double Fisher::p_value(void) const
  {
    if (calculate_p_exact())
      return p_value_exact();
    return p_value_approximative();
  }


  double Fisher::p_right(void) const
  {
    if (!calculate_p_exact()) {
      if (oddsratio_<1.0)
        return 1.0-p_value_approximative()/2.0;
      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_);
    assert(a_ <= a_+b_ && b_ <= a_+b_);
    assert(a_ <= a_+c_ && c_ <= a_+c_);

    return cdf_hypergeometric_P(c_, c_+d_, a_+b_, a_+c_);
  }


  double Fisher::p_value_one_sided(void) const
  {
    return p_right();
  }


  double Fisher::p_value_approximative(void) const
  {
    return gsl_cdf_chisq_Q(Chi2(), 1.0);
  }


  double Fisher::p_value_exact(void) const
  {
    // The hypergeometric distribution is unimodal (only one peak) and the
    // 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_);
    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);

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


  double Fisher::yates(unsigned int o, double e) const
  {
    double x = std::abs(o - e) - 0.5;
    return x*x/e;
  }

}}} // of namespace statistics, yat, and theplu