libs/math/test/test_hypergeometric_dist.cpp

   1 // Copyright John Maddock 2008
   2 // Copyright Paul A. Bristow
   3 // Copyright Gautam Sewani
   4
   5 // Use, modification and distribution are subject to the
   6 // Boost Software License, Version 1.0.
   7 // (See accompanying file LICENSE_1_0.txt
   8 // or copy at http://www.boost.org/LICENSE_1_0.txt)
   9
  10 #define BOOST_MATH_OVERFLOW_ERROR_POLICY throw_on_error
  11 #include <boost/math/concepts/real_concept.hpp> // for real_concept
  12 #include <boost/math/distributions/hypergeometric.hpp>
  13
  14 #define BOOST_TEST_MAIN
  15 #include <boost/test/unit_test.hpp> // Boost.Test
  16 #include <boost/test/results_collector.hpp>
  17 #include <boost/test/unit_test.hpp>
  18 #include <boost/test/tools/floating_point_comparison.hpp>
  19
  20 #include <iostream>
  21    using std::cout;
  22    using std::endl;
  23    using std::setprecision;
  24
  25 #include <boost/array.hpp>
  26 #include "functor.hpp"
  27 #include "handle_test_result.hpp"
  28 #include "table_type.hpp"
  29
  30 #define BOOST_CHECK_EX(a) \
  31    {\
  32       unsigned int failures = boost::unit_test::results_collector.results( boost::unit_test::framework::current_test_case().p_id ).p_assertions_failed;\
  33       BOOST_CHECK(a); \
  34       if(failures != boost::unit_test::results_collector.results( boost::unit_test::framework::current_test_case().p_id ).p_assertions_failed)\
  35       {\
  36          std::cerr << "Failure was with data ";\
  37          std::cerr << std::setprecision(35); \
  38          std::cerr << "x = " << x << ", r = " << r << ", n = " << n\
  39          << ", N = " << N << ", p = " << cp << ", q = " << ccp << std::endl;\
  40       }\
  41    }
  42
  43 void expected_results()
  44 {
  45    //
  46    // Define the max and mean errors expected for
  47    // various compilers and platforms.
  48    //
  49    const char* largest_type;
  50 #ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
  51    if(boost::math::policies::digits<double, boost::math::policies::policy<> >() == boost::math::policies::digits<long double, boost::math::policies::policy<> >())
  52    {
  53       largest_type = "(long\\s+)?double|real_concept";
  54    }
  55    else
  56    {
  57       largest_type = "long double|real_concept";
  58    }
  59 #else
  60    largest_type = "(long\\s+)?double";
  61 #endif
  62 #ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
  63    if((boost::math::tools::digits<long double>() > boost::math::tools::digits<double>())
  64       && (boost::math::tools::digits<long double>() < 100))
  65    {
  66       //
  67       // Some split of errors from long double into double:
  68       //
  69       add_expected_result(
  70          ".*",                          // compiler
  71          ".*",                          // stdlib
  72          ".*",                          // platform
  73          "double",                      // test type(s)
  74          "Random.*",                    // test data group
  75          ".*", 1500, 1500);      // test function
  76       add_expected_result(
  77          ".*",                          // compiler
  78          ".*",                          // stdlib
  79          ".*",                          // platform
  80          "double",                      // test type(s)
  81          ".*",                    // test data group
  82          ".*", 10, 10);      // test function
  83    }
  84 #endif
  85
  86    add_expected_result(
  87       ".*",                          // compiler
  88       ".*",                          // stdlib
  89       ".*",                          // platform
  90       "real_concept",                // test type(s)
  91       "Random.*",                    // test data group
  92       ".*", 250000000, 25000000);    // test function
  93    add_expected_result(
  94       ".*",                          // compiler
  95       ".*",                          // stdlib
  96       ".*",                          // platform
  97       largest_type,                  // test type(s)
  98       "Random.*",                    // test data group
  99       ".*", 10000000, 5000000);      // test function
 100    add_expected_result(
 101       ".*",                          // compiler
 102       ".*",                          // stdlib
 103       ".*",                          // platform
 104       largest_type,                  // test type(s)
 105       ".*",                          // test data group
 106       ".*", 50, 20);                 // test function
 107 }
 108
 109 template <class T>
 110 inline unsigned make_unsigned(T x)
 111 {
 112    return static_cast<unsigned>(x);
 113 }
 114 template<>
 115 inline unsigned make_unsigned(boost::math::concepts::real_concept x)
 116 {
 117    return static_cast<unsigned>(x.value());
 118 }
 119
 120 template <class T>
 121 T pdf_tester(T r, T n, T N, T x)
 122 {
 123    boost::math::hypergeometric_distribution<T> d(make_unsigned(r), make_unsigned(n), make_unsigned(N));
 124    return pdf(d, x);
 125 }
 126
 127 template <class T>
 128 T cdf_tester(T r, T n, T N, T x)
 129 {
 130    boost::math::hypergeometric_distribution<T> d(make_unsigned(r), make_unsigned(n), make_unsigned(N));
 131    return cdf(d, x);
 132 }
 133
 134 template <class T>
 135 T ccdf_tester(T r, T n, T N, T x)
 136 {
 137    boost::math::hypergeometric_distribution<T> d(make_unsigned(r), make_unsigned(n), make_unsigned(N));
 138    return cdf(complement(d, x));
 139 }
 140
 141 template <class Real, class T>
 142 void do_test_hypergeometric(const T& data, const char* type_name, const char* test_name)
 143 {
 144    // warning suppression:
 145    (void)data;
 146    (void)type_name;
 147    (void)test_name;
 148
 149 #if !defined(TEST_QUANT) || (TEST_QUANT == 0)
 150    typedef Real                   value_type;
 151
 152    typedef value_type (*pg)(value_type, value_type, value_type, value_type);
 153 #if defined(BOOST_MATH_NO_DEDUCED_FUNCTION_POINTERS)
 154    pg funcp = pdf_tester<value_type>;
 155 #else
 156    pg funcp = pdf_tester;
 157 #endif
 158
 159    boost::math::tools::test_result<value_type> result;
 160
 161    std::cout << "Testing " << test_name << " with type " << type_name
 162       << "\n~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\n";
 163
 164    //
 165    // test hypergeometric against data:
 166    //
 167    result = boost::math::tools::test_hetero<Real>(
 168       data,
 169       bind_func<Real>(funcp, 0, 1, 2, 3),
 170       extract_result<Real>(4));
 171    handle_test_result(result, data[result.worst()], result.worst(), type_name, "hypergeometric PDF", test_name);
 172
 173 #if defined(BOOST_MATH_NO_DEDUCED_FUNCTION_POINTERS)
 174    funcp = cdf_tester<value_type>;
 175 #else
 176    funcp = cdf_tester;
 177 #endif
 178
 179    //
 180    // test hypergeometric against data:
 181    //
 182    result = boost::math::tools::test_hetero<Real>(
 183       data,
 184       bind_func<Real>(funcp, 0, 1, 2, 3),
 185       extract_result<Real>(5));
 186    handle_test_result(result, data[result.worst()], result.worst(), type_name, "hypergeometric CDF", test_name);
 187
 188 #if defined(BOOST_MATH_NO_DEDUCED_FUNCTION_POINTERS)
 189    funcp = ccdf_tester<value_type>;
 190 #else
 191    funcp = ccdf_tester;
 192 #endif
 193
 194    //
 195    // test hypergeometric against data:
 196    //
 197    result = boost::math::tools::test_hetero<Real>(
 198       data,
 199       bind_func<Real>(funcp, 0, 1, 2, 3),
 200       extract_result<Real>(6));
 201    handle_test_result(result, data[result.worst()], result.worst(), type_name, "hypergeometric CDF complement", test_name);
 202    std::cout << std::endl;
 203 #endif
 204 }
 205
 206 template <class Real, class T>
 207 void do_test_hypergeometric_quantile(const T& data, const char* type_name, const char* test_name)
 208 {
 209    typedef Real                   value_type;
 210
 211    std::cout << "Checking quantiles with " << test_name << " with type " << type_name
 212       << "\n~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\n";
 213
 214    if(boost::math::tools::digits<value_type>() > 50)
 215    {
 216       for(unsigned i = 0; i < data.size(); ++i)
 217       {
 218          using namespace boost::math::policies;
 219
 220          unsigned r = make_unsigned(data[i][0]);
 221          unsigned n = make_unsigned(data[i][1]);
 222          unsigned N = make_unsigned(data[i][2]);
 223          unsigned x = make_unsigned(data[i][3]);
 224          value_type cp = data[i][5];
 225          value_type ccp = data[i][6];
 226          //
 227          // A bit of warning suppression:
 228          //
 229          (void)x;
 230          (void)n;
 231          (void)r;
 232          (void)N;
 233          (void)cp;
 234          (void)ccp;
 235
 236 #if !defined(TEST_QUANT) || (TEST_QUANT == 1)
 237          boost::math::hypergeometric_distribution<value_type,
 238             policy<discrete_quantile<integer_round_up> > > du(r, n, N);
 239
 240          if((cp < 0.9) && (cp > boost::math::tools::min_value<value_type>()))
 241          {
 242             BOOST_CHECK_EX(quantile(du, cp) >= x);
 243          }
 244          if((ccp < 0.9) && (ccp > boost::math::tools::min_value<value_type>()))
 245          {
 246             BOOST_CHECK_EX(quantile(complement(du, ccp)) >= x);
 247          }
 248 #endif
 249 #if !defined(TEST_QUANT) || (TEST_QUANT == 2)
 250          boost::math::hypergeometric_distribution<value_type,
 251             policy<discrete_quantile<integer_round_down> > > dl(r, n, N);
 252          if((cp < 0.9) && (cp > boost::math::tools::min_value<value_type>()))
 253          {
 254             BOOST_CHECK_EX(quantile(dl, cp) <= x);
 255          }
 256          if((ccp < 0.9) && (ccp > boost::math::tools::min_value<value_type>()))
 257          {
 258             BOOST_CHECK_EX(quantile(complement(dl, ccp)) <= x);
 259          }
 260 #endif
 261 #if !defined(TEST_QUANT) || (TEST_QUANT == 3)
 262          boost::math::hypergeometric_distribution<value_type,
 263             policy<discrete_quantile<integer_round_nearest> > > dn(r, n, N);
 264
 265          if((cp < 0.9) && (cp > boost::math::tools::min_value<value_type>()))
 266          {
 267             BOOST_CHECK_EX(quantile(dn, cp) == x);
 268          }
 269          if((ccp < 0.9) && (ccp > boost::math::tools::min_value<value_type>()))
 270          {
 271             BOOST_CHECK_EX(quantile(complement(dn, ccp)) == x);
 272          }
 273 #endif
 274 #if !defined(TEST_QUANT) || (TEST_QUANT == 4)
 275          boost::math::hypergeometric_distribution<value_type,
 276             policy<discrete_quantile<integer_round_outwards> > > dou(r, n, N);
 277
 278          if((cp < 0.9) && (cp > boost::math::tools::min_value<value_type>()))
 279          {
 280             if(cp < 0.5)
 281             {
 282                BOOST_CHECK_EX(quantile(dou, cp) <= x);
 283             }
 284             else
 285             {
 286                BOOST_CHECK_EX(quantile(dou, cp) >= x);
 287             }
 288          }
 289          if((ccp < 0.9) && (ccp > boost::math::tools::min_value<value_type>()))
 290          {
 291             if(ccp < 0.5)
 292             {
 293                BOOST_CHECK_EX(quantile(complement(dou, ccp)) >= x);
 294             }
 295             else
 296             {
 297                BOOST_CHECK_EX(quantile(complement(dou, ccp)) <= x);
 298             }
 299          }
 300 #endif
 301 #if !defined(TEST_QUANT) || (TEST_QUANT == 5)
 302          boost::math::hypergeometric_distribution<value_type,
 303             policy<discrete_quantile<integer_round_inwards> > > di(r, n, N);
 304
 305          if((cp < 0.9) && (cp > boost::math::tools::min_value<value_type>()))
 306          {
 307             if(cp < 0.5)
 308             {
 309                BOOST_CHECK_EX(quantile(di, cp) >= x);
 310             }
 311             else
 312             {
 313                BOOST_CHECK_EX(quantile(di, cp) <= x);
 314             }
 315          }
 316          if((ccp < 0.9) && (ccp > boost::math::tools::min_value<value_type>()))
 317          {
 318             if(ccp < 0.5)
 319             {
 320                BOOST_CHECK_EX(quantile(complement(di, ccp)) <= x);
 321             }
 322             else
 323             {
 324                BOOST_CHECK_EX(quantile(complement(di, ccp)) >= x);
 325             }
 326          }
 327 #endif
 328       }
 329    }
 330 }
 331
 332
 333 template <class RealType>
 334 void test_spot(unsigned x, unsigned n, unsigned r, unsigned N,
 335                RealType p, RealType cp, RealType ccp, RealType tol)
 336 {
 337    //
 338    // A bit of warning suppression:
 339    //
 340    (void)x;
 341    (void)n;
 342    (void)r;
 343    (void)N;
 344    (void)p;
 345    (void)cp;
 346    (void)ccp;
 347    (void)tol;
 348 #if !defined(TEST_QUANT) || (TEST_QUANT == 0)
 349    boost::math::hypergeometric_distribution<RealType> d(r, n, N);
 350
 351    std::pair<unsigned, unsigned> extent = range(d);
 352    // CDF's:
 353    BOOST_CHECK_CLOSE(pdf(d, x), p, tol);
 354    BOOST_CHECK_CLOSE(cdf(d, x), cp, tol);
 355    BOOST_CHECK_CLOSE(cdf(complement(d, x)), ccp, tol);
 356    // Again with real-value arguments:
 357    BOOST_CHECK_CLOSE(pdf(d, static_cast<RealType>(x)), p, tol);
 358    BOOST_CHECK_CLOSE(cdf(d, static_cast<RealType>(x)), cp, tol);
 359    BOOST_CHECK_CLOSE(cdf(complement(d, static_cast<RealType>(x))), ccp, tol);
 360    //
 361    // Quantiles, don't bother checking these for type float
 362    // as there's not enough precision in the p and q values
 363    // to get back to where we started:
 364    //
 365    if(boost::math::tools::digits<RealType>() > 50)
 366    {
 367       using namespace boost::math::policies;
 368
 369       boost::math::hypergeometric_distribution<RealType,
 370          policy<discrete_quantile<integer_round_up> > > du(r, n, N);
 371
 372       BOOST_CHECK_EX(quantile(du, cp) >= x);
 373       BOOST_CHECK_EX(quantile(complement(du, ccp)) >= x);
 374
 375       boost::math::hypergeometric_distribution<RealType,
 376          policy<discrete_quantile<integer_round_down> > > dl(r, n, N);
 377       BOOST_CHECK_EX(quantile(dl, cp) <= x);
 378       BOOST_CHECK_EX(quantile(complement(dl, ccp)) <= x);
 379
 380       boost::math::hypergeometric_distribution<RealType,
 381          policy<discrete_quantile<integer_round_nearest> > > dn(r, n, N);
 382
 383       BOOST_CHECK_EX(quantile(dn, cp) == x);
 384       BOOST_CHECK_EX(quantile(complement(dn, ccp)) == x);
 385    }
 386    //
 387    // Error checking of out of bounds arguments:
 388    //
 389    BOOST_MATH_CHECK_THROW(pdf(d, extent.second + 1), std::domain_error);
 390    BOOST_MATH_CHECK_THROW(cdf(d, extent.second + 1), std::domain_error);
 391    BOOST_MATH_CHECK_THROW(cdf(complement(d, extent.second + 1)), std::domain_error);
 392    if(extent.first > 0)
 393    {
 394       BOOST_MATH_CHECK_THROW(pdf(d, extent.first - 1), std::domain_error);
 395       BOOST_MATH_CHECK_THROW(cdf(d, extent.first - 1), std::domain_error);
 396       BOOST_MATH_CHECK_THROW(cdf(complement(d, extent.first - 1)), std::domain_error);
 397    }
 398    BOOST_MATH_CHECK_THROW(quantile(d, 1.1f), std::domain_error);
 399    BOOST_MATH_CHECK_THROW(quantile(complement(d, 1.1f)), std::domain_error);
 400    BOOST_MATH_CHECK_THROW(quantile(d, -0.001f), std::domain_error);
 401    BOOST_MATH_CHECK_THROW(quantile(complement(d, -0.001f)), std::domain_error);
 402    //
 403    // Checking of extreme values:
 404    //
 405    BOOST_CHECK_EQUAL(quantile(d, 0), extent.first);
 406    BOOST_CHECK_EQUAL(quantile(d, 1), extent.second);
 407    BOOST_CHECK_EQUAL(quantile(complement(d, 0)), extent.second);
 408    BOOST_CHECK_EQUAL(quantile(complement(d, 1)), extent.first);
 409    BOOST_CHECK_EQUAL(cdf(d, extent.second), 1);
 410    BOOST_CHECK_EQUAL(cdf(complement(d, extent.second)), 0);
 411 #endif
 412 }
 413
 414 template <class RealType>
 415 void test_spots(RealType /*T*/, const char* type_name)
 416 {
 417    // Basic sanity checks.
 418    // 50 eps as a percentage, up to a maximum of double precision
 419    // Test data taken from Mathematica 6
 420 #define T RealType
 421 #include "hypergeometric_test_data.ipp"
 422    do_test_hypergeometric<T>(hypergeometric_test_data, type_name, "Mathematica data");
 423
 424 #include "hypergeometric_dist_data2.ipp"
 425    if(boost::is_floating_point<RealType>::value)
 426    {
 427       //
 428       // Don't test this for real_concept: it's too slow!!!
 429       //
 430       do_test_hypergeometric<T>(hypergeometric_dist_data2, type_name, "Random large data");
 431    }
 432
 433    do_test_hypergeometric_quantile<T>(hypergeometric_test_data, type_name, "Mathematica data");
 434    if(boost::is_floating_point<RealType>::value)
 435    {
 436       //
 437       // Don't test this for real_concept: it's too slow!!!
 438       //
 439       do_test_hypergeometric_quantile<T>(hypergeometric_dist_data2, type_name, "Random large data");
 440    }
 441
 442    RealType tolerance = (std::max)(
 443       static_cast<RealType>(2e-16L),    // limit of test data
 444       boost::math::tools::epsilon<RealType>());
 445    cout<<"Absolute tolerance:"<<tolerance<<endl;
 446
 447    tolerance *= 50 * 100; // 50eps as a persentage
 448    cout << "Tolerance for type " << typeid(RealType).name()  << " is " << tolerance << " %" << endl;
 449
 450    //
 451    // These sanity check values were calculated using the online
 452    // calculator at http://stattrek.com/Tables/Hypergeometric.aspx
 453    // It's assumed that the test values are accurate to no more than
 454    // double precision.
 455    //
 456    test_spot(20, 200, 50, 500, static_cast<T>(0.120748236361163), static_cast<T>(0.563532430195156), static_cast<T>(1 - 0.563532430195156), tolerance);
 457    test_spot(53, 452, 64, 500, static_cast<T>(0.0184749573044286), static_cast<T>(0.0299118078796907), static_cast<T>(1 - 0.0299118078796907), tolerance);
 458    test_spot(32, 1287, 128, 5000, static_cast<T>(0.0807012167418264), static_cast<T>(0.469768774237742), static_cast<T>(1 - 0.469768774237742), tolerance);
 459    test_spot(1, 13, 4, 26, static_cast<T>(0.248695652173913), static_cast<T>(0.296521739130435), static_cast<T>(1 - 0.296521739130435), tolerance);
 460    test_spot(2, 13, 4, 26, static_cast<T>(0.40695652173913), static_cast<T>(0.703478260869565), static_cast<T>(1 - 0.703478260869565), tolerance);
 461    test_spot(3, 13, 4, 26, static_cast<T>(0.248695652173913), static_cast<T>(0.952173913043478), static_cast<T>(1 - 0.952173913043478), tolerance);
 462    test_spot(40, 70, 89, 170, static_cast<T>(0.0721901023798991), static_cast<T>(0.885447799131944), static_cast<T>(1 - 0.885447799131944), tolerance);
 463
 464    boost::math::hypergeometric_distribution<RealType> d(50, 200, 500);
 465    BOOST_CHECK_EQUAL(range(d).first, 0u);
 466    BOOST_CHECK_EQUAL(range(d).second, 50u);
 467    BOOST_CHECK_CLOSE(mean(d), static_cast<RealType>(20), tolerance);
 468    BOOST_CHECK_CLOSE(mode(d), static_cast<RealType>(20), tolerance);
 469    BOOST_CHECK_CLOSE(variance(d), static_cast<RealType>(10.821643286573146292585170340681L), tolerance);
 470    BOOST_CHECK_CLOSE(skewness(d), static_cast<RealType>(0.048833071022952084732902910189366L), tolerance);
 471    BOOST_CHECK_CLOSE(kurtosis_excess(d), static_cast<RealType>(2.5155486690782804816404001878293L), tolerance);
 472    BOOST_CHECK_CLOSE(kurtosis(d), kurtosis_excess(d) + 3, tolerance);
 473    BOOST_CHECK_EQUAL(quantile(d, 0.5f), median(d));
 474
 475    BOOST_MATH_CHECK_THROW(d = boost::math::hypergeometric_distribution<RealType>(501, 40, 500), std::domain_error);
 476    BOOST_MATH_CHECK_THROW(d = boost::math::hypergeometric_distribution<RealType>(40, 501, 500), std::domain_error);
 477 }
 478
 479
 480 BOOST_AUTO_TEST_CASE( test_main )
 481 {
 482    expected_results();
 483    // Basic sanity-check spot values.
 484    // (Parameter value, arbitrarily zero, only communicates the floating point type).
 485    test_spots(0.0F, "float"); // Test float. OK at decdigits = 0 tolerance = 0.0001 %
 486    test_spots(0.0, "double"); // Test double. OK at decdigits 7, tolerance = 1e07 %
 487 #ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
 488    test_spots(0.0L, "long double"); // Test long double.
 489 #ifndef BOOST_MATH_NO_REAL_CONCEPT_TESTS
 490    test_spots(boost::math::concepts::real_concept(0), "real_concept"); // Test real_concept.
 491 #endif
 492 #else
 493    std::cout << "<note>The long double tests have been disabled on this platform "
 494       "either because the long double overloads of the usual math functions are "
 495       "not available at all, or because they are too inaccurate for these tests "
 496       "to pass.</note>" << std::endl;
 497 #endif
 498
 499
 500 } // BOOST_AUTO_TEST_CASE( test_main )
 501