source: trunk/test/fisher.cc @ 3743

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

/*
  Copyright (C) 2008, 2009, 2010, 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 "Suite.h"

#include "yat/statistics/Fisher.h"

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

using namespace theplu::yat;
void test_p_value(test::Suite&);
void test_p_value_approximative(test::Suite&);
void test_p_value_exact(test::Suite&);
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[])
{
  test::Suite suite(argc, argv);

  statistics::Fisher f;
  if (!suite.equal(f.oddsratio(1,4,5,1), 0.05)) {
    suite.add(false);
    suite.err() << "oddsratio failed\n";
  }
  double a, b, c, d;
  f.expected(a,b,c,d);
  if (!suite.equal(a, 30.0/11.0)) {
    suite.add(false);
    suite.err() << "expected a failed\n";
  }
  if (!suite.equal(b, 25.0/11.0)) {
    suite.add(false);
    suite.err() << "expected b failed\n";
  }
  if (!suite.equal(c, 36.0/11.0)) {
    suite.add(false);
    suite.err() << "expected c failed\n";
  }
  if (!suite.equal(d, 30.0/11.0)) {
    suite.add(false);
    suite.err() << "expected d failed\n";
  }
  if (!suite.equal(f.Chi2(), 4.4122222222222222222222)) {
    suite.add(false);
    suite.err() << "Chi2 failed\n";
  }
  if (!suite.equal(f.minimum_size(),10)) {
    suite.add(false);
    suite.err() << "minimum_size failed\n";
  }
  test_p_value(suite);
  test_large_numbers(suite);
  test_yates(suite);
  test_cdf_hypergeometric_Q(suite);
  test_ticket908(suite);
  return suite.return_value();
}


void test_large_numbers(test::Suite& suite)
{
  // skip test if unsigned int is 16 bit
  if ((UINT_MAX >> 16) == 0) {
    suite.out() << "skipping test_large_numbers\n";
    return;
  }

  statistics::Fisher f;
  double oddsratio = f.oddsratio(1166,63326825-1166,1095,66074759-1095);
  if (oddsratio<0.5 || oddsratio>2) {
    suite.err() << "oddsratio: " << oddsratio << "\n";
    suite.err() << "expected ~ 1\n";
    suite.add(false);
  }
  suite.add(suite.equal_fix(f.p_value(), 0.0123, 0.0001));
  f.p_left();
  f.p_right();
}


void test_p_value(test::Suite& suite)
{
  test_p_value_exact(suite);
  test_p_value_approximative(suite);
}


void test_p_value_approximative(test::Suite& suite)
{
  suite.err() << "testing p_value_approximative\n";
  statistics::Fisher f;
  f.minimum_size() = 0;
  f.oddsratio(10,20,20,50);
  // oddsratio 1.25 > 1 so p = 2*p_right
  suite.add(suite.equal(f.p_value(), 2*f.p_right()));
  f.oddsratio(10,20,10,20);
  suite.add(suite.equal(f.p_value(), 1.0));
  suite.add(suite.equal(f.p_right(), 0.5));
  suite.add(suite.equal(f.p_left(), 0.5));
}


void test_p_value_exact(test::Suite& suite)
{
  suite.err() << "test p_value_exact\n";
  statistics::Fisher f;
  f.minimum_size() = 1000;
  f.oddsratio(10,20,10,20);
  suite.add(suite.equal(f.p_value(), 1.0, 20));

  f.oddsratio(10, 20, 20, 200);
  suite.add(suite.equal_fix(f.p_value(), 0.000811906062767622,1e-16));
  suite.add(suite.equal_fix(f.p_right(), 0.000811906062767622,1e-16));

  // testing symmetry
  statistics::Fisher f2;
  f2.minimum_size() = 1000;
  f2.oddsratio(20, 200, 10, 20);
  suite.add(suite.equal(f2.p_value(), f.p_value()));
  suite.add(suite.equal(f2.p_left(), f.p_right()));

  f.oddsratio(1, 1, 1, 2);
  suite.add(suite.equal(f.p_value(), 1.0, 2));
  suite.add(suite.equal(f.p_right(), 0.7));
  suite.add(suite.equal(f.p_left(), 0.9));

  f.oddsratio(1, 1, 2, 1);
  suite.add(suite.equal(f.p_value(), 1.0, 2));
  suite.add(suite.equal(f.p_right(), 0.9));
  suite.add(suite.equal(f.p_left(), 0.7));
}


void test_yates(test::Suite& suite)
{
  statistics::Fisher f;
  statistics::Fisher f2(true);

  f.oddsratio(5,10,10,10);
  f2.oddsratio(5,10,10,10);
  f.minimum_size() = 0;
  f2.minimum_size() = 0;
  double p = f.p_value();
  double p_yates = f2.p_value();
  if (p==p_yates) {
    suite.add(false);
    suite.err() << "p-value: " << p << "\n";
    suite.err() << "p-value yates's corrected: " << p_yates << "\n";
  }
}


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";
  }
}