Changeset 3743 for trunk/test/fisher.cc


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/test/fisher.cc

    r3455 r3743  
    22
    33/*
    4   Copyright (C) 2008, 2009, 2010, 2012, 2013, 2014, 2015 Peter Johansson
     4  Copyright (C) 2008, 2009, 2010, 2012, 2013, 2014, 2015, 2018 Peter Johansson
    55
    66  This file is part of the yat library, http://dev.thep.lu.se/yat
     
    2626#include "yat/statistics/Fisher.h"
    2727
     28#include <gsl/gsl_cdf.h>
     29#include <gsl/gsl_randist.h>
    2830#include <climits>
    2931
     
    3436void test_large_numbers(test::Suite&);
    3537void test_yates(test::Suite&);
     38void test_cdf_hypergeometric_Q(test::Suite&);
     39void test_ticket908(test::Suite&);
    3640
    3741int main(int argc, char* argv[])
     
    7377  test_large_numbers(suite);
    7478  test_yates(suite);
     79  test_cdf_hypergeometric_Q(suite);
     80  test_ticket908(suite);
    7581  return suite.return_value();
    7682}
     
    168174  }
    169175}
     176
     177
     178void test_cdf_hypergeometric_Q(test::Suite& suite)
     179{
     180  suite.out() << "\ntest " << __func__ << "\n";
     181  unsigned int n1 = 16;
     182  unsigned int n2 = 3;
     183  unsigned int t = 16;
     184  for (unsigned int k=0; k<=t; ++k) {
     185    // P(X=k)
     186    double p = gsl_ran_hypergeometric_pdf(k, n1, n2, t);
     187    // P(X<=k)
     188    double P = gsl_cdf_hypergeometric_P(k, n1, n2, t);
     189    // P(X>k)
     190    double Q = gsl_cdf_hypergeometric_Q(k, n1, n2, t);
     191    suite.out() << k << " " << p << " " << P << " " << Q << "\n";
     192    double sum = 0;
     193    for (unsigned int x=0; x<=k; ++x)
     194      sum += gsl_ran_hypergeometric_pdf(x, n1, n2, t);
     195    if (!suite.equal(sum, P, 20)) {
     196      suite.add(false);
     197      suite.err() << "error: gsl_cdf_hypergeometric_P for k=" << k << "\n";
     198    }
     199    if (!suite.equal(Q, 1-P, 1000)) {
     200      suite.add(false);
     201      suite.err() << "failed: Q = 1-P for k=" << k << "\n";
     202    }
     203  }
     204}
     205
     206
     207// test bug #908 reported in http://dev.thep.lu.se/yat/ticket/908
     208void test_ticket908(test::Suite& suite)
     209{
     210  suite.out() << "===\ntesting ticket 908\n";
     211  statistics::Fisher f;
     212  int a1 = 13;
     213  int a2 = 3;
     214  int b1 = 3;
     215  int b2 = 0;
     216
     217  double correct_p = 0;
     218  unsigned int n1 = a1 + a2;
     219  unsigned int n2 = b1 + b2;
     220  unsigned int t = a1+b1;
     221  for (unsigned int k=0; k<=n1; ++k) {
     222    double P = gsl_ran_hypergeometric_pdf(k, n1, n2, t);
     223    if (P <=  gsl_ran_hypergeometric_pdf(a1, n1, n2, t)) {
     224      suite.out() << "counting " << k << " " << P << "\n";
     225      correct_p += P;
     226    }
     227  }
     228  statistics::Fisher fisher;
     229  std::vector<double> p;
     230  fisher.oddsratio(a1, a2, b1, b2);
     231  p.push_back(fisher.p_value());
     232  fisher.oddsratio(a2, a1, b2, b1);
     233  p.push_back(fisher.p_value());
     234  fisher.oddsratio(b1, b2, a1, a2);
     235  p.push_back(fisher.p_value());
     236  fisher.oddsratio(b2, b1, a2, a1);
     237  p.push_back(fisher.p_value());
     238  for (size_t i=0; i<p.size(); ++i) {
     239    if (!suite.add(suite.equal(correct_p, p[i], 20)))
     240      suite.err() << "error: incorrect p: " << p[i] << " for case "
     241                  << i << "; expected: " << correct_p << "\n";
     242  }
     243}
Note: See TracChangeset for help on using the changeset viewer.