ceph/src/boost/libs/math/test/test_binomial.cpp

   1 // test_binomial.cpp
   2
   3 // Copyright John Maddock 2006.
   4 // Copyright  Paul A. Bristow 2007.
   5
   6 // Use, modification and distribution are subject to the
   7 // Boost Software License, Version 1.0.
   8 // (See accompanying file LICENSE_1_0.txt
   9 // or copy at http://www.boost.org/LICENSE_1_0.txt)
  10
  11 // Basic sanity test for Binomial Cumulative Distribution Function.
  12
  13 #define BOOST_MATH_DISCRETE_QUANTILE_POLICY real
  14
  15 #if !defined(TEST_FLOAT) && !defined(TEST_DOUBLE) && !defined(TEST_LDOUBLE) && !defined(TEST_REAL_CONCEPT)
  16 #  define TEST_FLOAT
  17 #  define TEST_DOUBLE
  18 #  define TEST_LDOUBLE
  19 #  define TEST_REAL_CONCEPT
  20 #endif
  21
  22 #ifdef _MSC_VER
  23 #  pragma warning(disable: 4127) // conditional expression is constant.
  24 #  pragma warning(disable: 4100) // unreferenced formal parameter.
  25 // Seems an entirely spurious warning - formal parameter T IS used - get error if /* T */
  26 //#  pragma warning(disable: 4535) // calling _set_se_translator() requires /EHa (in Boost.test)
  27 // Enable C++ Exceptions Yes With SEH Exceptions (/EHa) prevents warning 4535.
  28 #endif
  29
  30 #include <boost/math/tools/test.hpp>
  31 #include <boost/math/concepts/real_concept.hpp> // for real_concept
  32 using ::boost::math::concepts::real_concept;
  33
  34 #include <boost/math/distributions/binomial.hpp> // for binomial_distribution
  35 using boost::math::binomial_distribution;
  36
  37 #define BOOST_TEST_MAIN
  38 #include <boost/test/unit_test.hpp> // for test_main
  39 #include <boost/test/tools/floating_point_comparison.hpp> // for BOOST_CHECK_CLOSE
  40 #include "table_type.hpp"
  41
  42 #include "test_out_of_range.hpp"
  43
  44 #include <iostream>
  45 using std::cout;
  46 using std::endl;
  47 #include <limits>
  48 using std::numeric_limits;
  49
  50 template <class RealType>
  51 void test_spot(
  52      RealType N,    // Number of trials
  53      RealType k,    // Number of successes
  54      RealType p,    // Probability of success
  55      RealType P,    // CDF
  56      RealType Q,    // Complement of CDF
  57      RealType tol)  // Test tolerance
  58 {
  59    boost::math::binomial_distribution<RealType> bn(N, p);
  60    BOOST_CHECK_CLOSE(
  61       cdf(bn, k), P, tol);
  62    if((P < 0.99) && (Q < 0.99))
  63    {
  64       //
  65       // We can only check this if P is not too close to 1,
  66       // so that we can guarantee Q is free of error:
  67       //
  68       BOOST_CHECK_CLOSE(
  69          cdf(complement(bn, k)), Q, tol);
  70       if(k != 0)
  71       {
  72          BOOST_CHECK_CLOSE(
  73             quantile(bn, P), k, tol);
  74       }
  75       else
  76       {
  77          // Just check quantile is very small:
  78          if((std::numeric_limits<RealType>::max_exponent <= std::numeric_limits<double>::max_exponent) && (boost::is_floating_point<RealType>::value))
  79          {
  80             // Limit where this is checked: if exponent range is very large we may
  81             // run out of iterations in our root finding algorithm.
  82             BOOST_CHECK(quantile(bn, P) < boost::math::tools::epsilon<RealType>() * 10);
  83          }
  84       }
  85       if(k != 0)
  86       {
  87          BOOST_CHECK_CLOSE(
  88             quantile(complement(bn, Q)), k, tol);
  89       }
  90       else
  91       {
  92          // Just check quantile is very small:
  93          if((std::numeric_limits<RealType>::max_exponent <= std::numeric_limits<double>::max_exponent) && (boost::is_floating_point<RealType>::value))
  94          {
  95             // Limit where this is checked: if exponent range is very large we may
  96             // run out of iterations in our root finding algorithm.
  97             BOOST_CHECK(quantile(complement(bn, Q)) < boost::math::tools::epsilon<RealType>() * 10);
  98          }
  99       }
 100       if(k > 0)
 101       {
 102          // estimate success ratio:
 103          // Note lower bound uses a different formual internally
 104          // from upper bound, have to adjust things to prevent
 105          // fencepost errors:
 106          BOOST_CHECK_CLOSE(
 107             binomial_distribution<RealType>::find_lower_bound_on_p(
 108                N, k+1, Q),
 109             p, tol);
 110          BOOST_CHECK_CLOSE(
 111             binomial_distribution<RealType>::find_upper_bound_on_p(
 112                N, k, P),
 113             p, tol);
 114
 115          if(Q < P)
 116          {
 117             // Default method (Clopper Pearson)
 118             BOOST_CHECK(
 119                binomial_distribution<RealType>::find_lower_bound_on_p(
 120                   N, k, Q)
 121                   <=
 122                binomial_distribution<RealType>::find_upper_bound_on_p(
 123                   N, k, Q)
 124                   );
 125             BOOST_CHECK((
 126                binomial_distribution<RealType>::find_lower_bound_on_p(
 127                   N, k, Q)
 128                   <= k/N) && (k/N <=
 129                binomial_distribution<RealType>::find_upper_bound_on_p(
 130                   N, k, Q))
 131                   );
 132             // Bayes Method (Jeffreys Prior)
 133             BOOST_CHECK(
 134                binomial_distribution<RealType>::find_lower_bound_on_p(
 135                N, k, Q, binomial_distribution<RealType>::jeffreys_prior_interval)
 136                   <=
 137                binomial_distribution<RealType>::find_upper_bound_on_p(
 138                   N, k, Q, binomial_distribution<RealType>::jeffreys_prior_interval)
 139                   );
 140             BOOST_CHECK((
 141                binomial_distribution<RealType>::find_lower_bound_on_p(
 142                   N, k, Q, binomial_distribution<RealType>::jeffreys_prior_interval)
 143                   <= k/N) && (k/N <=
 144                binomial_distribution<RealType>::find_upper_bound_on_p(
 145                   N, k, Q, binomial_distribution<RealType>::jeffreys_prior_interval))
 146                   );
 147          }
 148          else
 149          {
 150             // Default method (Clopper Pearson)
 151             BOOST_CHECK(
 152                binomial_distribution<RealType>::find_lower_bound_on_p(
 153                   N, k, P)
 154                   <=
 155                binomial_distribution<RealType>::find_upper_bound_on_p(
 156                   N, k, P)
 157                   );
 158             BOOST_CHECK(
 159                (binomial_distribution<RealType>::find_lower_bound_on_p(
 160                   N, k, P)
 161                   <= k / N) && (k/N <=
 162                binomial_distribution<RealType>::find_upper_bound_on_p(
 163                   N, k, P))
 164                   );
 165             // Bayes Method (Jeffreys Prior)
 166             BOOST_CHECK(
 167                binomial_distribution<RealType>::find_lower_bound_on_p(
 168                   N, k, P, binomial_distribution<RealType>::jeffreys_prior_interval)
 169                   <=
 170                binomial_distribution<RealType>::find_upper_bound_on_p(
 171                   N, k, P, binomial_distribution<RealType>::jeffreys_prior_interval)
 172                   );
 173             BOOST_CHECK(
 174                (binomial_distribution<RealType>::find_lower_bound_on_p(
 175                   N, k, P, binomial_distribution<RealType>::jeffreys_prior_interval)
 176                   <= k / N) && (k/N <=
 177                binomial_distribution<RealType>::find_upper_bound_on_p(
 178                   N, k, P, binomial_distribution<RealType>::jeffreys_prior_interval))
 179                   );
 180          }
 181       }
 182       //
 183       // estimate sample size:
 184       //
 185       BOOST_CHECK_CLOSE(
 186          binomial_distribution<RealType>::find_minimum_number_of_trials(
 187             k, p, P),
 188          N, tol);
 189       BOOST_CHECK_CLOSE(
 190          binomial_distribution<RealType>::find_maximum_number_of_trials(
 191             k, p, Q),
 192          N, tol);
 193    }
 194
 195    // Double check consistency of CDF and PDF by computing
 196    // the finite sum:
 197    RealType sum = 0;
 198    for(unsigned i = 0; i <= k; ++i)
 199       sum += pdf(bn, RealType(i));
 200    BOOST_CHECK_CLOSE(
 201       sum, P, tol);
 202    // And complement as well:
 203    sum = 0;
 204    for(RealType i = N; i > k; i -= 1)
 205       sum += pdf(bn, i);
 206    if(P < 0.99)
 207    {
 208       BOOST_CHECK_CLOSE(
 209          sum, Q, tol);
 210    }
 211    else
 212    {
 213       // Not enough information content in P for Q to be meaningful
 214       RealType tol = (std::max)(2 * Q, boost::math::tools::epsilon<RealType>());
 215       BOOST_CHECK(sum < tol);
 216    }
 217 }
 218
 219 template <class RealType> // Any floating-point type RealType.
 220 void test_spots(RealType T)
 221 {
 222   // Basic sanity checks, test data is to double precision only
 223   // so set tolerance to 100eps expressed as a persent, or
 224   // 100eps of type double expressed as a persent, whichever
 225   // is the larger.
 226
 227   RealType tolerance = (std::max)
 228       (boost::math::tools::epsilon<RealType>(),
 229       static_cast<RealType>(std::numeric_limits<double>::epsilon()));
 230   tolerance *= 100 * 1000;
 231   RealType tol2 = boost::math::tools::epsilon<RealType>() * 5 * 100;  // 5 eps as a persent
 232
 233   cout << "Tolerance for type " << typeid(T).name()  << " is " << tolerance << " %" << endl;
 234
 235
 236   // Sources of spot test values:
 237
 238   // MathCAD defines pbinom(k, n, p)
 239   // returns pr(X ,=k) when random variable X has the binomial distribution with parameters n and p.
 240   // 0 <= k ,= n
 241   // 0 <= p <= 1
 242   // P = pbinom(30, 500, 0.05) = 0.869147702104609
 243
 244   using boost::math::binomial_distribution;
 245   using  ::boost::math::cdf;
 246   using  ::boost::math::pdf;
 247
 248 #if !defined(TEST_ROUNDING) || (TEST_ROUNDING == 0)
 249   // Test binomial using cdf spot values from MathCAD.
 250   // These test quantiles and complements as well.
 251   test_spot(
 252      static_cast<RealType>(500),                     // Sample size, N
 253      static_cast<RealType>(30),                      // Number of successes, k
 254      static_cast<RealType>(0.05),                    // Probability of success, p
 255      static_cast<RealType>(0.869147702104609),       // Probability of result (CDF), P
 256      static_cast<RealType>(1 - 0.869147702104609),   // Q = 1 - P
 257      tolerance);
 258
 259   test_spot(
 260      static_cast<RealType>(500),                     // Sample size, N
 261      static_cast<RealType>(250),                     // Number of successes, k
 262      static_cast<RealType>(0.05),                    // Probability of success, p
 263      static_cast<RealType>(1),                       // Probability of result (CDF), P
 264      static_cast<RealType>(0),   // Q = 1 - P
 265      tolerance);
 266
 267   test_spot(
 268      static_cast<RealType>(500),                     // Sample size, N
 269      static_cast<RealType>(470),                     // Number of successes, k
 270      static_cast<RealType>(0.95),                    // Probability of success, p
 271      static_cast<RealType>(0.176470742656766),       // Probability of result (CDF), P
 272      static_cast<RealType>(1 - 0.176470742656766),   // Q = 1 - P
 273      tolerance * 10);                                // Note higher tolerance on this test!
 274
 275   test_spot(
 276      static_cast<RealType>(500),                       // Sample size, N
 277      static_cast<RealType>(400),                       // Number of successes, k
 278      static_cast<RealType>(0.05),                      // Probability of success, p
 279      static_cast<RealType>(1),                         // Probability of result (CDF), P
 280      static_cast<RealType>(0),                         // Q = 1 - P
 281      tolerance);
 282
 283   test_spot(
 284      static_cast<RealType>(500),                       // Sample size, N
 285      static_cast<RealType>(400),                       // Number of successes, k
 286      static_cast<RealType>(0.9),                       // Probability of success, p
 287      static_cast<RealType>(1.80180425681923E-11),      // Probability of result (CDF), P
 288      static_cast<RealType>(1 - 1.80180425681923E-11),  // Q = 1 - P
 289      tolerance);
 290
 291   test_spot(
 292      static_cast<RealType>(500),                       // Sample size, N
 293      static_cast<RealType>(5),                         // Number of successes, k
 294      static_cast<RealType>(0.05),                      // Probability of success, p
 295      static_cast<RealType>(9.181808267643E-7),         // Probability of result (CDF), P
 296      static_cast<RealType>(1 - 9.181808267643E-7),     // Q = 1 - P
 297      tolerance);
 298
 299   test_spot(
 300      static_cast<RealType>(2),                       // Sample size, N
 301      static_cast<RealType>(1),                       // Number of successes, k
 302      static_cast<RealType>(0.5),                     // Probability of success, p
 303      static_cast<RealType>(0.75),                    // Probability of result (CDF), P
 304      static_cast<RealType>(0.25),                    // Q = 1 - P
 305      tolerance);
 306
 307   test_spot(
 308      static_cast<RealType>(8),                       // Sample size, N
 309      static_cast<RealType>(3),                       // Number of successes, k
 310      static_cast<RealType>(0.25),                    // Probability of success, p
 311      static_cast<RealType>(0.8861846923828125),      // Probability of result (CDF), P
 312      static_cast<RealType>(1 - 0.8861846923828125),  // Q = 1 - P
 313      tolerance);
 314
 315   test_spot(
 316      static_cast<RealType>(8),                       // Sample size, N
 317      static_cast<RealType>(0),                       // Number of successes, k
 318      static_cast<RealType>(0.25),                    // Probability of success, p
 319      static_cast<RealType>(0.1001129150390625),      // Probability of result (CDF), P
 320      static_cast<RealType>(1 - 0.1001129150390625),  // Q = 1 - P
 321      tolerance);
 322
 323   test_spot(
 324      static_cast<RealType>(8),                       // Sample size, N
 325      static_cast<RealType>(1),                       // Number of successes, k
 326      static_cast<RealType>(0.25),                    // Probability of success, p
 327      static_cast<RealType>(0.36708068847656244),     // Probability of result (CDF), P
 328      static_cast<RealType>(1 - 0.36708068847656244), // Q = 1 - P
 329      tolerance);
 330
 331   test_spot(
 332      static_cast<RealType>(8),                       // Sample size, N
 333      static_cast<RealType>(4),                       // Number of successes, k
 334      static_cast<RealType>(0.25),                    // Probability of success, p
 335      static_cast<RealType>(0.9727020263671875),      // Probability of result (CDF), P
 336      static_cast<RealType>(1 - 0.9727020263671875),  // Q = 1 - P
 337      tolerance);
 338
 339   test_spot(
 340      static_cast<RealType>(8),                       // Sample size, N
 341      static_cast<RealType>(7),                       // Number of successes, k
 342      static_cast<RealType>(0.25),                    // Probability of success, p
 343      static_cast<RealType>(0.9999847412109375),      // Probability of result (CDF), P
 344      static_cast<RealType>(1 - 0.9999847412109375),  // Q = 1 - P
 345      tolerance);
 346
 347   // Tests on PDF follow:
 348   BOOST_CHECK_CLOSE(
 349      pdf(binomial_distribution<RealType>(static_cast<RealType>(20), static_cast<RealType>(0.75)),
 350      static_cast<RealType>(10)),  // k.
 351      static_cast<RealType>(0.00992227527967770583927631378173), // 0.00992227527967770583927631378173
 352      tolerance);
 353
 354   BOOST_CHECK_CLOSE(
 355     pdf(binomial_distribution<RealType>(static_cast<RealType>(20), static_cast<RealType>(0.5)),
 356     static_cast<RealType>(10)),  // k.
 357     static_cast<RealType>(0.17619705200195312500000000000000000000), // get k=10 0.049611376398388612 p = 0.25
 358     tolerance);
 359
 360   // Binomial pdf Test values from
 361   // http://www.adsciengineering.com/bpdcalc/index.php  for example
 362   // http://www.adsciengineering.com/bpdcalc/index.php?n=20&p=0.25&start=0&stop=20&Submit=Generate
 363   // Appears to use at least 80-bit long double for 32 decimal digits accuracy,
 364   // but loses accuracy of display if leading zeros?
 365   // (if trailings zero then are exact values?)
 366   // so useful for testing 64-bit double accuracy.
 367   // P = 0.25, n = 20, k = 0 to 20
 368
 369   //0   C(20,0) * 0.25^0 * 0.75^20   0.00317121193893399322405457496643
 370   //1   C(20,1) * 0.25^1 * 0.75^19   0.02114141292622662149369716644287
 371   //2   C(20,2) * 0.25^2 * 0.75^18   0.06694780759971763473004102706909
 372   //3   C(20,3) * 0.25^3 * 0.75^17   0.13389561519943526946008205413818
 373   //4   C(20,4) * 0.25^4 * 0.75^16   0.18968545486586663173511624336242
 374   //5   C(20,5) * 0.25^5 * 0.75^15   0.20233115185692440718412399291992
 375   //6   C(20,6) * 0.25^6 * 0.75^14   0.16860929321410367265343666076660
 376   //7   C(20,7) * 0.25^7 * 0.75^13   0.11240619547606911510229110717773
 377   //8   C(20,8) * 0.25^8 * 0.75^12   0.06088668921620410401374101638793
 378   //9   C(20,9) * 0.25^9 * 0.75^11   0.02706075076275737956166267395019
 379   //10   C(20,10) * 0.25^10 * 0.75^10   0.00992227527967770583927631378173
 380   //11   C(20,11) * 0.25^11 * 0.75^9   0.00300675008475081995129585266113
 381   //12   C(20,12) * 0.25^12 * 0.75^8   0.00075168752118770498782396316528
 382   //13   C(20,13) * 0.25^13 * 0.75^7   0.00015419231203850358724594116210
 383   //14   C(20,14) * 0.25^14 * 0.75^6   0.00002569871867308393120765686035
 384   //15   C(20,15) * 0.25^15 * 0.75^5   0.00000342649582307785749435424804
 385   //16   C(20,16) * 0.25^16 * 0.75^4   0.00000035692664823727682232856750
 386   //17   C(20,17) * 0.25^17 * 0.75^3   0.00000002799424692057073116302490
 387   //18   C(20,18) * 0.25^18 * 0.75^2   0.00000000155523594003170728683471
 388   //19   C(20,19) * 0.25^19 * 0.75^1   0.00000000005456968210637569427490
 389   //20   C(20,20) * 0.25^20 * 0.75^0   0.00000000000090949470177292823791
 390
 391
 392     BOOST_CHECK_CLOSE(
 393     pdf(binomial_distribution<RealType>(static_cast<RealType>(20), static_cast<RealType>(0.25)),
 394     static_cast<RealType>(10)),  // k.
 395     static_cast<RealType>(0.00992227527967770583927631378173), // k=10  p = 0.25
 396     tolerance);
 397
 398     BOOST_CHECK_CLOSE( // k = 0 use different formula - only exp so more accurate.
 399     pdf(binomial_distribution<RealType>(static_cast<RealType>(20), static_cast<RealType>(0.25)),
 400     static_cast<RealType>(0)),  // k.
 401     static_cast<RealType>(0.00317121193893399322405457496643), // k=0  p = 0.25
 402     tolerance);
 403
 404     BOOST_CHECK_CLOSE( // k = 20 use different formula - only exp so more accurate.
 405     pdf(binomial_distribution<RealType>(static_cast<RealType>(20), static_cast<RealType>(0.25)),
 406     static_cast<RealType>(20)),  // k == n.
 407     static_cast<RealType>(0.00000000000090949470177292823791), // k=20  p = 0.25
 408     tolerance);
 409
 410     BOOST_CHECK_CLOSE( // k = 1.
 411     pdf(binomial_distribution<RealType>(static_cast<RealType>(20), static_cast<RealType>(0.25)),
 412     static_cast<RealType>(1)),  // k.
 413     static_cast<RealType>(0.02114141292622662149369716644287), // k=1  p = 0.25
 414     tolerance);
 415
 416     // Some exact (probably) values.
 417     BOOST_CHECK_CLOSE(
 418     pdf(binomial_distribution<RealType>(static_cast<RealType>(8), static_cast<RealType>(0.25)),
 419     static_cast<RealType>(0)),  // k.
 420     static_cast<RealType>(0.10011291503906250000000000000000), // k=0  p = 0.25
 421     tolerance);
 422
 423     BOOST_CHECK_CLOSE( // k = 1.
 424     pdf(binomial_distribution<RealType>(static_cast<RealType>(8), static_cast<RealType>(0.25)),
 425     static_cast<RealType>(1)),  // k.
 426     static_cast<RealType>(0.26696777343750000000000000000000), // k=1  p = 0.25
 427     tolerance);
 428
 429     BOOST_CHECK_CLOSE( // k = 2.
 430     pdf(binomial_distribution<RealType>(static_cast<RealType>(8), static_cast<RealType>(0.25)),
 431     static_cast<RealType>(2)),  // k.
 432     static_cast<RealType>(0.31146240234375000000000000000000), // k=2  p = 0.25
 433     tolerance);
 434
 435     BOOST_CHECK_CLOSE( // k = 3.
 436     pdf(binomial_distribution<RealType>(static_cast<RealType>(8), static_cast<RealType>(0.25)),
 437     static_cast<RealType>(3)),  // k.
 438     static_cast<RealType>(0.20764160156250000000000000000000), // k=3  p = 0.25
 439     tolerance);
 440
 441     BOOST_CHECK_CLOSE( // k = 7.
 442     pdf(binomial_distribution<RealType>(static_cast<RealType>(8), static_cast<RealType>(0.25)),
 443     static_cast<RealType>(7)),  // k.
 444     static_cast<RealType>(0.00036621093750000000000000000000), // k=7  p = 0.25
 445     tolerance);
 446
 447     BOOST_CHECK_CLOSE( // k = 8.
 448     pdf(binomial_distribution<RealType>(static_cast<RealType>(8), static_cast<RealType>(0.25)),
 449     static_cast<RealType>(8)),  // k = n.
 450     static_cast<RealType>(0.00001525878906250000000000000000), // k=8  p = 0.25
 451     tolerance);
 452
 453     binomial_distribution<RealType> dist(static_cast<RealType>(8), static_cast<RealType>(0.25));
 454     RealType x = static_cast<RealType>(0.125);
 455     using namespace std; // ADL of std names.
 456     // mean:
 457     BOOST_CHECK_CLOSE(
 458        mean(dist)
 459        , static_cast<RealType>(8 * 0.25), tol2);
 460     // variance:
 461     BOOST_CHECK_CLOSE(
 462        variance(dist)
 463        , static_cast<RealType>(8 * 0.25 * 0.75), tol2);
 464     // std deviation:
 465     BOOST_CHECK_CLOSE(
 466        standard_deviation(dist)
 467        , static_cast<RealType>(sqrt(8 * 0.25L * 0.75L)), tol2);
 468     // hazard:
 469     BOOST_CHECK_CLOSE(
 470        hazard(dist, x)
 471        , pdf(dist, x) / cdf(complement(dist, x)), tol2);
 472     // cumulative hazard:
 473     BOOST_CHECK_CLOSE(
 474        chf(dist, x)
 475        , -log(cdf(complement(dist, x))), tol2);
 476     // coefficient_of_variation:
 477     BOOST_CHECK_CLOSE(
 478        coefficient_of_variation(dist)
 479        , standard_deviation(dist) / mean(dist), tol2);
 480     // mode:
 481     BOOST_CHECK_CLOSE(
 482        mode(dist)
 483        , static_cast<RealType>(std::floor(9 * 0.25)), tol2);
 484     // skewness:
 485     BOOST_CHECK_CLOSE(
 486        skewness(dist)
 487        , static_cast<RealType>(0.40824829046386301636621401245098L), (std::max)(tol2, static_cast<RealType>(5e-29))); // test data has 32 digits only.
 488     // kurtosis:
 489     BOOST_CHECK_CLOSE(
 490        kurtosis(dist)
 491        , static_cast<RealType>(2.916666666666666666666666666666666666L), tol2);
 492     // kurtosis excess:
 493     BOOST_CHECK_CLOSE(
 494        kurtosis_excess(dist)
 495        , static_cast<RealType>(-0.08333333333333333333333333333333333333L), tol2);
 496     // Check kurtosis_excess == kurtosis -3;
 497       BOOST_CHECK_EQUAL(kurtosis(dist), static_cast<RealType>(3) + kurtosis_excess(dist));
 498
 499     // special cases for PDF:
 500     BOOST_CHECK_EQUAL(
 501        pdf(
 502           binomial_distribution<RealType>(static_cast<RealType>(8), static_cast<RealType>(0)),
 503           static_cast<RealType>(0)), static_cast<RealType>(1)
 504        );
 505     BOOST_CHECK_EQUAL(
 506        pdf(
 507           binomial_distribution<RealType>(static_cast<RealType>(8), static_cast<RealType>(0)),
 508           static_cast<RealType>(0.0001)), static_cast<RealType>(0)
 509        );
 510     BOOST_CHECK_EQUAL(
 511        pdf(
 512           binomial_distribution<RealType>(static_cast<RealType>(8), static_cast<RealType>(1)),
 513           static_cast<RealType>(0.001)), static_cast<RealType>(0)
 514        );
 515     BOOST_CHECK_EQUAL(
 516        pdf(
 517           binomial_distribution<RealType>(static_cast<RealType>(8), static_cast<RealType>(1)),
 518           static_cast<RealType>(8)), static_cast<RealType>(1)
 519        );
 520     BOOST_CHECK_EQUAL(
 521        pdf(
 522           binomial_distribution<RealType>(static_cast<RealType>(0), static_cast<RealType>(0.25)),
 523           static_cast<RealType>(0)), static_cast<RealType>(1)
 524        );
 525     BOOST_MATH_CHECK_THROW(
 526        pdf(
 527           binomial_distribution<RealType>(static_cast<RealType>(-1), static_cast<RealType>(0.25)),
 528           static_cast<RealType>(0)), std::domain_error
 529        );
 530     BOOST_MATH_CHECK_THROW(
 531        pdf(
 532           binomial_distribution<RealType>(static_cast<RealType>(8), static_cast<RealType>(-0.25)),
 533           static_cast<RealType>(0)), std::domain_error
 534        );
 535     BOOST_MATH_CHECK_THROW(
 536        pdf(
 537           binomial_distribution<RealType>(static_cast<RealType>(8), static_cast<RealType>(1.25)),
 538           static_cast<RealType>(0)), std::domain_error
 539        );
 540     BOOST_MATH_CHECK_THROW(
 541        pdf(
 542           binomial_distribution<RealType>(static_cast<RealType>(8), static_cast<RealType>(0.25)),
 543           static_cast<RealType>(-1)), std::domain_error
 544        );
 545     BOOST_MATH_CHECK_THROW(
 546        pdf(
 547           binomial_distribution<RealType>(static_cast<RealType>(8), static_cast<RealType>(0.25)),
 548           static_cast<RealType>(9)), std::domain_error
 549        );
 550     BOOST_MATH_CHECK_THROW(
 551        cdf(
 552           binomial_distribution<RealType>(static_cast<RealType>(8), static_cast<RealType>(0.25)),
 553           static_cast<RealType>(-1)), std::domain_error
 554        );
 555     BOOST_MATH_CHECK_THROW(
 556        cdf(
 557           binomial_distribution<RealType>(static_cast<RealType>(8), static_cast<RealType>(0.25)),
 558           static_cast<RealType>(9)), std::domain_error
 559        );
 560     BOOST_MATH_CHECK_THROW(
 561        cdf(
 562           binomial_distribution<RealType>(static_cast<RealType>(8), static_cast<RealType>(-0.25)),
 563           static_cast<RealType>(0)), std::domain_error
 564        );
 565     BOOST_MATH_CHECK_THROW(
 566        cdf(
 567           binomial_distribution<RealType>(static_cast<RealType>(8), static_cast<RealType>(1.25)),
 568           static_cast<RealType>(0)), std::domain_error
 569        );
 570     BOOST_MATH_CHECK_THROW(
 571        quantile(
 572           binomial_distribution<RealType>(static_cast<RealType>(8), static_cast<RealType>(-0.25)),
 573           static_cast<RealType>(0)), std::domain_error
 574        );
 575     BOOST_MATH_CHECK_THROW(
 576        quantile(
 577           binomial_distribution<RealType>(static_cast<RealType>(8), static_cast<RealType>(1.25)),
 578           static_cast<RealType>(0)), std::domain_error
 579        );
 580
 581     BOOST_CHECK_EQUAL(
 582        quantile(
 583           binomial_distribution<RealType>(static_cast<RealType>(16), static_cast<RealType>(0.25)),
 584           static_cast<RealType>(0.01)), // Less than cdf == pdf(binomial_distribution<RealType>(16, 0.25), 0)
 585           static_cast<RealType>(0) // so expect zero as best approximation.
 586        );
 587
 588     BOOST_CHECK_EQUAL(
 589        cdf(
 590           binomial_distribution<RealType>(static_cast<RealType>(8), static_cast<RealType>(0.25)),
 591           static_cast<RealType>(8)), static_cast<RealType>(1)
 592        );
 593     BOOST_CHECK_EQUAL(
 594        cdf(
 595           binomial_distribution<RealType>(static_cast<RealType>(8), static_cast<RealType>(0)),
 596           static_cast<RealType>(7)), static_cast<RealType>(1)
 597        );
 598     BOOST_CHECK_EQUAL(
 599        cdf(
 600           binomial_distribution<RealType>(static_cast<RealType>(8), static_cast<RealType>(1)),
 601           static_cast<RealType>(7)), static_cast<RealType>(0)
 602        );
 603
 604 #endif
 605
 606   {
 607     // This is a visual sanity check that everything is OK:
 608     binomial_distribution<RealType> my8dist(8., 0.25); // Note: double values (matching the distribution definition) avoid the need for any casting.
 609     //cout << "mean(my8dist) = " << boost::math::mean(my8dist) << endl; // mean(my8dist) = 2
 610     //cout << "my8dist.trials() = " << my8dist.trials()  << endl; // my8dist.trials() = 8
 611     //cout << "my8dist.success_fraction() = " << my8dist.success_fraction()  << endl; // my8dist.success_fraction() = 0.25
 612     BOOST_CHECK_CLOSE(my8dist.trials(), static_cast<RealType>(8), tol2);
 613     BOOST_CHECK_CLOSE(my8dist.success_fraction(), static_cast<RealType>(0.25), tol2);
 614
 615    //{
 616    //   int n = static_cast<int>(boost::math::tools::real_cast<double>(my8dist.trials()));
 617    //   RealType sumcdf = 0.;
 618    //   for (int k = 0; k <= n; k++)
 619    //   {
 620    //     cout << k << ' ' << pdf(my8dist, static_cast<RealType>(k));
 621    //     sumcdf += pdf(my8dist, static_cast<RealType>(k));
 622    //     cout  << ' '  << sumcdf;
 623    //     cout << ' ' << cdf(my8dist, static_cast<RealType>(k));
 624    //     cout << ' ' << sumcdf - cdf(my8dist, static_cast<RealType>(k)) << endl;
 625    //   } // for k
 626    // }
 627     // n = 8, p =0.25
 628     //k         pdf              cdf
 629     //0 0.1001129150390625 0.1001129150390625
 630     //1 0.26696777343749994 0.36708068847656244
 631     //2 0.31146240234375017 0.67854309082031261
 632     //3 0.20764160156249989 0.8861846923828125
 633     //4 0.086517333984375 0.9727020263671875
 634     //5 0.023071289062499997 0.9957733154296875
 635     //6 0.0038452148437500009 0.9996185302734375
 636     //7 0.00036621093749999984 0.9999847412109375
 637     //8 1.52587890625e-005 1 1 0
 638   }
 639 #define T RealType
 640 #include "binomial_quantile.ipp"
 641
 642   for(unsigned i = 0; i < binomial_quantile_data.size(); ++i)
 643   {
 644      using namespace boost::math::policies;
 645      RealType tol = boost::math::tools::epsilon<RealType>() * 500;
 646      if(!boost::is_floating_point<RealType>::value)
 647         tol *= 10;  // no lanczos approximation implies less accuracy
 648      RealType x;
 649 #if !defined(TEST_ROUNDING) || (TEST_ROUNDING == 1)
 650      //
 651      // Check full real value first:
 652      //
 653      typedef policy<discrete_quantile<boost::math::policies::real> > P1;
 654      binomial_distribution<RealType, P1> p1(binomial_quantile_data[i][0], binomial_quantile_data[i][1]);
 655      x = quantile(p1, binomial_quantile_data[i][2]);
 656      BOOST_CHECK_CLOSE_FRACTION(x, (RealType)binomial_quantile_data[i][3], tol);
 657      x = quantile(complement(p1, (RealType)binomial_quantile_data[i][2]));
 658      BOOST_CHECK_CLOSE_FRACTION(x, (RealType)binomial_quantile_data[i][4], tol);
 659 #endif
 660 #if !defined(TEST_ROUNDING) || (TEST_ROUNDING == 2)
 661      //
 662      // Now with round down to integer:
 663      //
 664      typedef policy<discrete_quantile<integer_round_down> > P2;
 665      binomial_distribution<RealType, P2> p2(binomial_quantile_data[i][0], binomial_quantile_data[i][1]);
 666      x = quantile(p2, binomial_quantile_data[i][2]);
 667      BOOST_CHECK_EQUAL(x, (RealType)floor(binomial_quantile_data[i][3]));
 668      x = quantile(complement(p2, binomial_quantile_data[i][2]));
 669      BOOST_CHECK_EQUAL(x, (RealType)floor(binomial_quantile_data[i][4]));
 670 #endif
 671 #if !defined(TEST_ROUNDING) || (TEST_ROUNDING == 3)
 672      //
 673      // Now with round up to integer:
 674      //
 675      typedef policy<discrete_quantile<integer_round_up> > P3;
 676      binomial_distribution<RealType, P3> p3(binomial_quantile_data[i][0], binomial_quantile_data[i][1]);
 677      x = quantile(p3, binomial_quantile_data[i][2]);
 678      BOOST_CHECK_EQUAL(x, (RealType)ceil(binomial_quantile_data[i][3]));
 679      x = quantile(complement(p3, binomial_quantile_data[i][2]));
 680      BOOST_CHECK_EQUAL(x, (RealType)ceil(binomial_quantile_data[i][4]));
 681 #endif
 682 #if !defined(TEST_ROUNDING) || (TEST_ROUNDING == 4)
 683      //
 684      // Now with round to integer "outside":
 685      //
 686      typedef policy<discrete_quantile<integer_round_outwards> > P4;
 687      binomial_distribution<RealType, P4> p4(binomial_quantile_data[i][0], binomial_quantile_data[i][1]);
 688      x = quantile(p4, binomial_quantile_data[i][2]);
 689      BOOST_CHECK_EQUAL(x, (RealType)(binomial_quantile_data[i][2] < 0.5f ? floor(binomial_quantile_data[i][3]) : ceil(binomial_quantile_data[i][3])));
 690      x = quantile(complement(p4, binomial_quantile_data[i][2]));
 691      BOOST_CHECK_EQUAL(x, (RealType)(binomial_quantile_data[i][2] < 0.5f ? ceil(binomial_quantile_data[i][4]) : floor(binomial_quantile_data[i][4])));
 692 #endif
 693 #if !defined(TEST_ROUNDING) || (TEST_ROUNDING == 5)
 694      //
 695      // Now with round to integer "inside":
 696      //
 697      typedef policy<discrete_quantile<integer_round_inwards> > P5;
 698      binomial_distribution<RealType, P5> p5(binomial_quantile_data[i][0], binomial_quantile_data[i][1]);
 699      x = quantile(p5, binomial_quantile_data[i][2]);
 700      BOOST_CHECK_EQUAL(x, (RealType)(binomial_quantile_data[i][2] < 0.5f ? ceil(binomial_quantile_data[i][3]) : floor(binomial_quantile_data[i][3])));
 701      x = quantile(complement(p5, binomial_quantile_data[i][2]));
 702      BOOST_CHECK_EQUAL(x, (RealType)(binomial_quantile_data[i][2] < 0.5f ? floor(binomial_quantile_data[i][4]) : ceil(binomial_quantile_data[i][4])));
 703 #endif
 704 #if !defined(TEST_ROUNDING) || (TEST_ROUNDING == 6)
 705      //
 706      // Now with round to nearest integer:
 707      //
 708      typedef policy<discrete_quantile<integer_round_nearest> > P6;
 709      binomial_distribution<RealType, P6> p6(binomial_quantile_data[i][0], binomial_quantile_data[i][1]);
 710      x = quantile(p6, binomial_quantile_data[i][2]);
 711      BOOST_CHECK_EQUAL(x, (RealType)(floor(binomial_quantile_data[i][3] + 0.5f)));
 712      x = quantile(complement(p6, binomial_quantile_data[i][2]));
 713      BOOST_CHECK_EQUAL(x, (RealType)(floor(binomial_quantile_data[i][4] + 0.5f)));
 714 #endif
 715   }
 716
 717    check_out_of_range<boost::math::binomial_distribution<RealType> >(1, 1); // (All) valid constructor parameter values.
 718
 719
 720 } // template <class RealType>void test_spots(RealType)
 721
 722 BOOST_AUTO_TEST_CASE( test_main )
 723 {
 724    BOOST_MATH_CONTROL_FP;
 725    // Check that can generate binomial distribution using one convenience methods:
 726    binomial_distribution<> mybn2(1., 0.5); // Using default RealType double.
 727   // but that
 728    // boost::math::binomial mybn1(1., 0.5); // Using typedef fails
 729   // error C2039: 'binomial' : is not a member of 'boost::math'
 730
 731   // Basic sanity-check spot values.
 732
 733   // (Parameter value, arbitrarily zero, only communicates the floating point type).
 734 #ifdef TEST_FLOAT
 735   test_spots(0.0F); // Test float.
 736 #endif
 737 #ifdef TEST_DOUBLE
 738   test_spots(0.0); // Test double.
 739 #endif
 740 #ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
 741 #ifdef TEST_LDOUBLE
 742   test_spots(0.0L); // Test long double.
 743 #endif
 744 #if !defined(BOOST_MATH_NO_REAL_CONCEPT_TESTS)
 745 #ifdef TEST_REAL_CONCEPT
 746   test_spots(boost::math::concepts::real_concept(0.)); // Test real concept.
 747 #endif
 748 #endif
 749 #else
 750    std::cout << "<note>The long double tests have been disabled on this platform "
 751       "either because the long double overloads of the usual math functions are "
 752       "not available at all, or because they are too inaccurate for these tests "
 753       "to pass.</note>" << std::endl;
 754 #endif
 755
 756 } // BOOST_AUTO_TEST_CASE( test_main )
 757
 758 /*
 759
 760 Output is:
 761
 762   Description: Autorun "J:\Cpp\MathToolkit\test\Math_test\Debug\test_binomial.exe"
 763   Running 1 test case...
 764   Tolerance for type float is 0.0119209 %
 765   Tolerance for type double is 2.22045e-011 %
 766   Tolerance for type long double is 2.22045e-011 %
 767   Tolerance for type class boost::math::concepts::real_concept is 2.22045e-011 %
 768
 769   *** No errors detected
 770
 771 ========== Build: 1 succeeded, 0 failed, 0 up-to-date, 0 skipped ==========
 772
 773
 774 */