123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700701702703704705706707708709710711712713714715716717718719720721722723724725726727728729730731732733734735736737738739740741742743744745746747748749750751752753754755756757758759760761762763764765766767768769770771772773774 |
- // test_binomial.cpp
- // Copyright John Maddock 2006.
- // Copyright Paul A. Bristow 2007.
- // Use, modification and distribution are subject to the
- // Boost Software License, Version 1.0.
- // (See accompanying file LICENSE_1_0.txt
- // or copy at http://www.boost.org/LICENSE_1_0.txt)
- // Basic sanity test for Binomial Cumulative Distribution Function.
- #define BOOST_MATH_DISCRETE_QUANTILE_POLICY real
- #if !defined(TEST_FLOAT) && !defined(TEST_DOUBLE) && !defined(TEST_LDOUBLE) && !defined(TEST_REAL_CONCEPT)
- # define TEST_FLOAT
- # define TEST_DOUBLE
- # define TEST_LDOUBLE
- # define TEST_REAL_CONCEPT
- #endif
- #ifdef _MSC_VER
- # pragma warning(disable: 4127) // conditional expression is constant.
- # pragma warning(disable: 4100) // unreferenced formal parameter.
- // Seems an entirely spurious warning - formal parameter T IS used - get error if /* T */
- //# pragma warning(disable: 4535) // calling _set_se_translator() requires /EHa (in Boost.test)
- // Enable C++ Exceptions Yes With SEH Exceptions (/EHa) prevents warning 4535.
- #endif
- #include <boost/math/tools/test.hpp>
- #include <boost/math/concepts/real_concept.hpp> // for real_concept
- using ::boost::math::concepts::real_concept;
- #include <boost/math/distributions/binomial.hpp> // for binomial_distribution
- using boost::math::binomial_distribution;
- #define BOOST_TEST_MAIN
- #include <boost/test/unit_test.hpp> // for test_main
- #include <boost/test/tools/floating_point_comparison.hpp> // for BOOST_CHECK_CLOSE
- #include "table_type.hpp"
- #include "test_out_of_range.hpp"
- #include <iostream>
- using std::cout;
- using std::endl;
- #include <limits>
- using std::numeric_limits;
- template <class RealType>
- void test_spot(
- RealType N, // Number of trials
- RealType k, // Number of successes
- RealType p, // Probability of success
- RealType P, // CDF
- RealType Q, // Complement of CDF
- RealType tol) // Test tolerance
- {
- boost::math::binomial_distribution<RealType> bn(N, p);
- BOOST_CHECK_CLOSE(
- cdf(bn, k), P, tol);
- if((P < 0.99) && (Q < 0.99))
- {
- //
- // We can only check this if P is not too close to 1,
- // so that we can guarantee Q is free of error:
- //
- BOOST_CHECK_CLOSE(
- cdf(complement(bn, k)), Q, tol);
- if(k != 0)
- {
- BOOST_CHECK_CLOSE(
- quantile(bn, P), k, tol);
- }
- else
- {
- // Just check quantile is very small:
- if((std::numeric_limits<RealType>::max_exponent <= std::numeric_limits<double>::max_exponent) && (boost::is_floating_point<RealType>::value))
- {
- // Limit where this is checked: if exponent range is very large we may
- // run out of iterations in our root finding algorithm.
- BOOST_CHECK(quantile(bn, P) < boost::math::tools::epsilon<RealType>() * 10);
- }
- }
- if(k != 0)
- {
- BOOST_CHECK_CLOSE(
- quantile(complement(bn, Q)), k, tol);
- }
- else
- {
- // Just check quantile is very small:
- if((std::numeric_limits<RealType>::max_exponent <= std::numeric_limits<double>::max_exponent) && (boost::is_floating_point<RealType>::value))
- {
- // Limit where this is checked: if exponent range is very large we may
- // run out of iterations in our root finding algorithm.
- BOOST_CHECK(quantile(complement(bn, Q)) < boost::math::tools::epsilon<RealType>() * 10);
- }
- }
- if(k > 0)
- {
- // estimate success ratio:
- // Note lower bound uses a different formual internally
- // from upper bound, have to adjust things to prevent
- // fencepost errors:
- BOOST_CHECK_CLOSE(
- binomial_distribution<RealType>::find_lower_bound_on_p(
- N, k+1, Q),
- p, tol);
- BOOST_CHECK_CLOSE(
- binomial_distribution<RealType>::find_upper_bound_on_p(
- N, k, P),
- p, tol);
- if(Q < P)
- {
- // Default method (Clopper Pearson)
- BOOST_CHECK(
- binomial_distribution<RealType>::find_lower_bound_on_p(
- N, k, Q)
- <=
- binomial_distribution<RealType>::find_upper_bound_on_p(
- N, k, Q)
- );
- BOOST_CHECK((
- binomial_distribution<RealType>::find_lower_bound_on_p(
- N, k, Q)
- <= k/N) && (k/N <=
- binomial_distribution<RealType>::find_upper_bound_on_p(
- N, k, Q))
- );
- // Bayes Method (Jeffreys Prior)
- BOOST_CHECK(
- binomial_distribution<RealType>::find_lower_bound_on_p(
- N, k, Q, binomial_distribution<RealType>::jeffreys_prior_interval)
- <=
- binomial_distribution<RealType>::find_upper_bound_on_p(
- N, k, Q, binomial_distribution<RealType>::jeffreys_prior_interval)
- );
- BOOST_CHECK((
- binomial_distribution<RealType>::find_lower_bound_on_p(
- N, k, Q, binomial_distribution<RealType>::jeffreys_prior_interval)
- <= k/N) && (k/N <=
- binomial_distribution<RealType>::find_upper_bound_on_p(
- N, k, Q, binomial_distribution<RealType>::jeffreys_prior_interval))
- );
- }
- else
- {
- // Default method (Clopper Pearson)
- BOOST_CHECK(
- binomial_distribution<RealType>::find_lower_bound_on_p(
- N, k, P)
- <=
- binomial_distribution<RealType>::find_upper_bound_on_p(
- N, k, P)
- );
- BOOST_CHECK(
- (binomial_distribution<RealType>::find_lower_bound_on_p(
- N, k, P)
- <= k / N) && (k/N <=
- binomial_distribution<RealType>::find_upper_bound_on_p(
- N, k, P))
- );
- // Bayes Method (Jeffreys Prior)
- BOOST_CHECK(
- binomial_distribution<RealType>::find_lower_bound_on_p(
- N, k, P, binomial_distribution<RealType>::jeffreys_prior_interval)
- <=
- binomial_distribution<RealType>::find_upper_bound_on_p(
- N, k, P, binomial_distribution<RealType>::jeffreys_prior_interval)
- );
- BOOST_CHECK(
- (binomial_distribution<RealType>::find_lower_bound_on_p(
- N, k, P, binomial_distribution<RealType>::jeffreys_prior_interval)
- <= k / N) && (k/N <=
- binomial_distribution<RealType>::find_upper_bound_on_p(
- N, k, P, binomial_distribution<RealType>::jeffreys_prior_interval))
- );
- }
- }
- //
- // estimate sample size:
- //
- BOOST_CHECK_CLOSE(
- binomial_distribution<RealType>::find_minimum_number_of_trials(
- k, p, P),
- N, tol);
- BOOST_CHECK_CLOSE(
- binomial_distribution<RealType>::find_maximum_number_of_trials(
- k, p, Q),
- N, tol);
- }
- // Double check consistency of CDF and PDF by computing
- // the finite sum:
- RealType sum = 0;
- for(unsigned i = 0; i <= k; ++i)
- sum += pdf(bn, RealType(i));
- BOOST_CHECK_CLOSE(
- sum, P, tol);
- // And complement as well:
- sum = 0;
- for(RealType i = N; i > k; i -= 1)
- sum += pdf(bn, i);
- if(P < 0.99)
- {
- BOOST_CHECK_CLOSE(
- sum, Q, tol);
- }
- else
- {
- // Not enough information content in P for Q to be meaningful
- RealType tol = (std::max)(2 * Q, boost::math::tools::epsilon<RealType>());
- BOOST_CHECK(sum < tol);
- }
- }
- template <class RealType> // Any floating-point type RealType.
- void test_spots(RealType T)
- {
- // Basic sanity checks, test data is to double precision only
- // so set tolerance to 100eps expressed as a persent, or
- // 100eps of type double expressed as a persent, whichever
- // is the larger.
- RealType tolerance = (std::max)
- (boost::math::tools::epsilon<RealType>(),
- static_cast<RealType>(std::numeric_limits<double>::epsilon()));
- tolerance *= 100 * 1000;
- RealType tol2 = boost::math::tools::epsilon<RealType>() * 5 * 100; // 5 eps as a persent
- cout << "Tolerance for type " << typeid(T).name() << " is " << tolerance << " %" << endl;
- // Sources of spot test values:
- // MathCAD defines pbinom(k, n, p)
- // returns pr(X ,=k) when random variable X has the binomial distribution with parameters n and p.
- // 0 <= k ,= n
- // 0 <= p <= 1
- // P = pbinom(30, 500, 0.05) = 0.869147702104609
- using boost::math::binomial_distribution;
- using ::boost::math::cdf;
- using ::boost::math::pdf;
- #if !defined(TEST_ROUNDING) || (TEST_ROUNDING == 0)
- // Test binomial using cdf spot values from MathCAD.
- // These test quantiles and complements as well.
- test_spot(
- static_cast<RealType>(500), // Sample size, N
- static_cast<RealType>(30), // Number of successes, k
- static_cast<RealType>(0.05), // Probability of success, p
- static_cast<RealType>(0.869147702104609), // Probability of result (CDF), P
- static_cast<RealType>(1 - 0.869147702104609), // Q = 1 - P
- tolerance);
- test_spot(
- static_cast<RealType>(500), // Sample size, N
- static_cast<RealType>(250), // Number of successes, k
- static_cast<RealType>(0.05), // Probability of success, p
- static_cast<RealType>(1), // Probability of result (CDF), P
- static_cast<RealType>(0), // Q = 1 - P
- tolerance);
- test_spot(
- static_cast<RealType>(500), // Sample size, N
- static_cast<RealType>(470), // Number of successes, k
- static_cast<RealType>(0.95), // Probability of success, p
- static_cast<RealType>(0.176470742656766), // Probability of result (CDF), P
- static_cast<RealType>(1 - 0.176470742656766), // Q = 1 - P
- tolerance * 10); // Note higher tolerance on this test!
- test_spot(
- static_cast<RealType>(500), // Sample size, N
- static_cast<RealType>(400), // Number of successes, k
- static_cast<RealType>(0.05), // Probability of success, p
- static_cast<RealType>(1), // Probability of result (CDF), P
- static_cast<RealType>(0), // Q = 1 - P
- tolerance);
- test_spot(
- static_cast<RealType>(500), // Sample size, N
- static_cast<RealType>(400), // Number of successes, k
- static_cast<RealType>(0.9), // Probability of success, p
- static_cast<RealType>(1.80180425681923E-11), // Probability of result (CDF), P
- static_cast<RealType>(1 - 1.80180425681923E-11), // Q = 1 - P
- tolerance);
- test_spot(
- static_cast<RealType>(500), // Sample size, N
- static_cast<RealType>(5), // Number of successes, k
- static_cast<RealType>(0.05), // Probability of success, p
- static_cast<RealType>(9.181808267643E-7), // Probability of result (CDF), P
- static_cast<RealType>(1 - 9.181808267643E-7), // Q = 1 - P
- tolerance);
- test_spot(
- static_cast<RealType>(2), // Sample size, N
- static_cast<RealType>(1), // Number of successes, k
- static_cast<RealType>(0.5), // Probability of success, p
- static_cast<RealType>(0.75), // Probability of result (CDF), P
- static_cast<RealType>(0.25), // Q = 1 - P
- tolerance);
- test_spot(
- static_cast<RealType>(8), // Sample size, N
- static_cast<RealType>(3), // Number of successes, k
- static_cast<RealType>(0.25), // Probability of success, p
- static_cast<RealType>(0.8861846923828125), // Probability of result (CDF), P
- static_cast<RealType>(1 - 0.8861846923828125), // Q = 1 - P
- tolerance);
- test_spot(
- static_cast<RealType>(8), // Sample size, N
- static_cast<RealType>(0), // Number of successes, k
- static_cast<RealType>(0.25), // Probability of success, p
- static_cast<RealType>(0.1001129150390625), // Probability of result (CDF), P
- static_cast<RealType>(1 - 0.1001129150390625), // Q = 1 - P
- tolerance);
- test_spot(
- static_cast<RealType>(8), // Sample size, N
- static_cast<RealType>(1), // Number of successes, k
- static_cast<RealType>(0.25), // Probability of success, p
- static_cast<RealType>(0.36708068847656244), // Probability of result (CDF), P
- static_cast<RealType>(1 - 0.36708068847656244), // Q = 1 - P
- tolerance);
- test_spot(
- static_cast<RealType>(8), // Sample size, N
- static_cast<RealType>(4), // Number of successes, k
- static_cast<RealType>(0.25), // Probability of success, p
- static_cast<RealType>(0.9727020263671875), // Probability of result (CDF), P
- static_cast<RealType>(1 - 0.9727020263671875), // Q = 1 - P
- tolerance);
- test_spot(
- static_cast<RealType>(8), // Sample size, N
- static_cast<RealType>(7), // Number of successes, k
- static_cast<RealType>(0.25), // Probability of success, p
- static_cast<RealType>(0.9999847412109375), // Probability of result (CDF), P
- static_cast<RealType>(1 - 0.9999847412109375), // Q = 1 - P
- tolerance);
- // Tests on PDF follow:
- BOOST_CHECK_CLOSE(
- pdf(binomial_distribution<RealType>(static_cast<RealType>(20), static_cast<RealType>(0.75)),
- static_cast<RealType>(10)), // k.
- static_cast<RealType>(0.00992227527967770583927631378173), // 0.00992227527967770583927631378173
- tolerance);
- BOOST_CHECK_CLOSE(
- pdf(binomial_distribution<RealType>(static_cast<RealType>(20), static_cast<RealType>(0.5)),
- static_cast<RealType>(10)), // k.
- static_cast<RealType>(0.17619705200195312500000000000000000000), // get k=10 0.049611376398388612 p = 0.25
- tolerance);
- // Binomial pdf Test values from
- // http://www.adsciengineering.com/bpdcalc/index.php for example
- // http://www.adsciengineering.com/bpdcalc/index.php?n=20&p=0.25&start=0&stop=20&Submit=Generate
- // Appears to use at least 80-bit long double for 32 decimal digits accuracy,
- // but loses accuracy of display if leading zeros?
- // (if trailings zero then are exact values?)
- // so useful for testing 64-bit double accuracy.
- // P = 0.25, n = 20, k = 0 to 20
- //0 C(20,0) * 0.25^0 * 0.75^20 0.00317121193893399322405457496643
- //1 C(20,1) * 0.25^1 * 0.75^19 0.02114141292622662149369716644287
- //2 C(20,2) * 0.25^2 * 0.75^18 0.06694780759971763473004102706909
- //3 C(20,3) * 0.25^3 * 0.75^17 0.13389561519943526946008205413818
- //4 C(20,4) * 0.25^4 * 0.75^16 0.18968545486586663173511624336242
- //5 C(20,5) * 0.25^5 * 0.75^15 0.20233115185692440718412399291992
- //6 C(20,6) * 0.25^6 * 0.75^14 0.16860929321410367265343666076660
- //7 C(20,7) * 0.25^7 * 0.75^13 0.11240619547606911510229110717773
- //8 C(20,8) * 0.25^8 * 0.75^12 0.06088668921620410401374101638793
- //9 C(20,9) * 0.25^9 * 0.75^11 0.02706075076275737956166267395019
- //10 C(20,10) * 0.25^10 * 0.75^10 0.00992227527967770583927631378173
- //11 C(20,11) * 0.25^11 * 0.75^9 0.00300675008475081995129585266113
- //12 C(20,12) * 0.25^12 * 0.75^8 0.00075168752118770498782396316528
- //13 C(20,13) * 0.25^13 * 0.75^7 0.00015419231203850358724594116210
- //14 C(20,14) * 0.25^14 * 0.75^6 0.00002569871867308393120765686035
- //15 C(20,15) * 0.25^15 * 0.75^5 0.00000342649582307785749435424804
- //16 C(20,16) * 0.25^16 * 0.75^4 0.00000035692664823727682232856750
- //17 C(20,17) * 0.25^17 * 0.75^3 0.00000002799424692057073116302490
- //18 C(20,18) * 0.25^18 * 0.75^2 0.00000000155523594003170728683471
- //19 C(20,19) * 0.25^19 * 0.75^1 0.00000000005456968210637569427490
- //20 C(20,20) * 0.25^20 * 0.75^0 0.00000000000090949470177292823791
- BOOST_CHECK_CLOSE(
- pdf(binomial_distribution<RealType>(static_cast<RealType>(20), static_cast<RealType>(0.25)),
- static_cast<RealType>(10)), // k.
- static_cast<RealType>(0.00992227527967770583927631378173), // k=10 p = 0.25
- tolerance);
- BOOST_CHECK_CLOSE( // k = 0 use different formula - only exp so more accurate.
- pdf(binomial_distribution<RealType>(static_cast<RealType>(20), static_cast<RealType>(0.25)),
- static_cast<RealType>(0)), // k.
- static_cast<RealType>(0.00317121193893399322405457496643), // k=0 p = 0.25
- tolerance);
- BOOST_CHECK_CLOSE( // k = 20 use different formula - only exp so more accurate.
- pdf(binomial_distribution<RealType>(static_cast<RealType>(20), static_cast<RealType>(0.25)),
- static_cast<RealType>(20)), // k == n.
- static_cast<RealType>(0.00000000000090949470177292823791), // k=20 p = 0.25
- tolerance);
- BOOST_CHECK_CLOSE( // k = 1.
- pdf(binomial_distribution<RealType>(static_cast<RealType>(20), static_cast<RealType>(0.25)),
- static_cast<RealType>(1)), // k.
- static_cast<RealType>(0.02114141292622662149369716644287), // k=1 p = 0.25
- tolerance);
- // Some exact (probably) values.
- BOOST_CHECK_CLOSE(
- pdf(binomial_distribution<RealType>(static_cast<RealType>(8), static_cast<RealType>(0.25)),
- static_cast<RealType>(0)), // k.
- static_cast<RealType>(0.10011291503906250000000000000000), // k=0 p = 0.25
- tolerance);
- BOOST_CHECK_CLOSE( // k = 1.
- pdf(binomial_distribution<RealType>(static_cast<RealType>(8), static_cast<RealType>(0.25)),
- static_cast<RealType>(1)), // k.
- static_cast<RealType>(0.26696777343750000000000000000000), // k=1 p = 0.25
- tolerance);
- BOOST_CHECK_CLOSE( // k = 2.
- pdf(binomial_distribution<RealType>(static_cast<RealType>(8), static_cast<RealType>(0.25)),
- static_cast<RealType>(2)), // k.
- static_cast<RealType>(0.31146240234375000000000000000000), // k=2 p = 0.25
- tolerance);
- BOOST_CHECK_CLOSE( // k = 3.
- pdf(binomial_distribution<RealType>(static_cast<RealType>(8), static_cast<RealType>(0.25)),
- static_cast<RealType>(3)), // k.
- static_cast<RealType>(0.20764160156250000000000000000000), // k=3 p = 0.25
- tolerance);
- BOOST_CHECK_CLOSE( // k = 7.
- pdf(binomial_distribution<RealType>(static_cast<RealType>(8), static_cast<RealType>(0.25)),
- static_cast<RealType>(7)), // k.
- static_cast<RealType>(0.00036621093750000000000000000000), // k=7 p = 0.25
- tolerance);
- BOOST_CHECK_CLOSE( // k = 8.
- pdf(binomial_distribution<RealType>(static_cast<RealType>(8), static_cast<RealType>(0.25)),
- static_cast<RealType>(8)), // k = n.
- static_cast<RealType>(0.00001525878906250000000000000000), // k=8 p = 0.25
- tolerance);
- binomial_distribution<RealType> dist(static_cast<RealType>(8), static_cast<RealType>(0.25));
- RealType x = static_cast<RealType>(0.125);
- using namespace std; // ADL of std names.
- // mean:
- BOOST_CHECK_CLOSE(
- mean(dist)
- , static_cast<RealType>(8 * 0.25), tol2);
- // variance:
- BOOST_CHECK_CLOSE(
- variance(dist)
- , static_cast<RealType>(8 * 0.25 * 0.75), tol2);
- // std deviation:
- BOOST_CHECK_CLOSE(
- standard_deviation(dist)
- , static_cast<RealType>(sqrt(8 * 0.25L * 0.75L)), tol2);
- // hazard:
- BOOST_CHECK_CLOSE(
- hazard(dist, x)
- , pdf(dist, x) / cdf(complement(dist, x)), tol2);
- // cumulative hazard:
- BOOST_CHECK_CLOSE(
- chf(dist, x)
- , -log(cdf(complement(dist, x))), tol2);
- // coefficient_of_variation:
- BOOST_CHECK_CLOSE(
- coefficient_of_variation(dist)
- , standard_deviation(dist) / mean(dist), tol2);
- // mode:
- BOOST_CHECK_CLOSE(
- mode(dist)
- , static_cast<RealType>(std::floor(9 * 0.25)), tol2);
- // skewness:
- BOOST_CHECK_CLOSE(
- skewness(dist)
- , static_cast<RealType>(0.40824829046386301636621401245098L), (std::max)(tol2, static_cast<RealType>(5e-29))); // test data has 32 digits only.
- // kurtosis:
- BOOST_CHECK_CLOSE(
- kurtosis(dist)
- , static_cast<RealType>(2.916666666666666666666666666666666666L), tol2);
- // kurtosis excess:
- BOOST_CHECK_CLOSE(
- kurtosis_excess(dist)
- , static_cast<RealType>(-0.08333333333333333333333333333333333333L), tol2);
- // Check kurtosis_excess == kurtosis -3;
- BOOST_CHECK_EQUAL(kurtosis(dist), static_cast<RealType>(3) + kurtosis_excess(dist));
- // special cases for PDF:
- BOOST_CHECK_EQUAL(
- pdf(
- binomial_distribution<RealType>(static_cast<RealType>(8), static_cast<RealType>(0)),
- static_cast<RealType>(0)), static_cast<RealType>(1)
- );
- BOOST_CHECK_EQUAL(
- pdf(
- binomial_distribution<RealType>(static_cast<RealType>(8), static_cast<RealType>(0)),
- static_cast<RealType>(0.0001)), static_cast<RealType>(0)
- );
- BOOST_CHECK_EQUAL(
- pdf(
- binomial_distribution<RealType>(static_cast<RealType>(8), static_cast<RealType>(1)),
- static_cast<RealType>(0.001)), static_cast<RealType>(0)
- );
- BOOST_CHECK_EQUAL(
- pdf(
- binomial_distribution<RealType>(static_cast<RealType>(8), static_cast<RealType>(1)),
- static_cast<RealType>(8)), static_cast<RealType>(1)
- );
- BOOST_CHECK_EQUAL(
- pdf(
- binomial_distribution<RealType>(static_cast<RealType>(0), static_cast<RealType>(0.25)),
- static_cast<RealType>(0)), static_cast<RealType>(1)
- );
- BOOST_MATH_CHECK_THROW(
- pdf(
- binomial_distribution<RealType>(static_cast<RealType>(-1), static_cast<RealType>(0.25)),
- static_cast<RealType>(0)), std::domain_error
- );
- BOOST_MATH_CHECK_THROW(
- pdf(
- binomial_distribution<RealType>(static_cast<RealType>(8), static_cast<RealType>(-0.25)),
- static_cast<RealType>(0)), std::domain_error
- );
- BOOST_MATH_CHECK_THROW(
- pdf(
- binomial_distribution<RealType>(static_cast<RealType>(8), static_cast<RealType>(1.25)),
- static_cast<RealType>(0)), std::domain_error
- );
- BOOST_MATH_CHECK_THROW(
- pdf(
- binomial_distribution<RealType>(static_cast<RealType>(8), static_cast<RealType>(0.25)),
- static_cast<RealType>(-1)), std::domain_error
- );
- BOOST_MATH_CHECK_THROW(
- pdf(
- binomial_distribution<RealType>(static_cast<RealType>(8), static_cast<RealType>(0.25)),
- static_cast<RealType>(9)), std::domain_error
- );
- BOOST_MATH_CHECK_THROW(
- cdf(
- binomial_distribution<RealType>(static_cast<RealType>(8), static_cast<RealType>(0.25)),
- static_cast<RealType>(-1)), std::domain_error
- );
- BOOST_MATH_CHECK_THROW(
- cdf(
- binomial_distribution<RealType>(static_cast<RealType>(8), static_cast<RealType>(0.25)),
- static_cast<RealType>(9)), std::domain_error
- );
- BOOST_MATH_CHECK_THROW(
- cdf(
- binomial_distribution<RealType>(static_cast<RealType>(8), static_cast<RealType>(-0.25)),
- static_cast<RealType>(0)), std::domain_error
- );
- BOOST_MATH_CHECK_THROW(
- cdf(
- binomial_distribution<RealType>(static_cast<RealType>(8), static_cast<RealType>(1.25)),
- static_cast<RealType>(0)), std::domain_error
- );
- BOOST_MATH_CHECK_THROW(
- quantile(
- binomial_distribution<RealType>(static_cast<RealType>(8), static_cast<RealType>(-0.25)),
- static_cast<RealType>(0)), std::domain_error
- );
- BOOST_MATH_CHECK_THROW(
- quantile(
- binomial_distribution<RealType>(static_cast<RealType>(8), static_cast<RealType>(1.25)),
- static_cast<RealType>(0)), std::domain_error
- );
- BOOST_CHECK_EQUAL(
- quantile(
- binomial_distribution<RealType>(static_cast<RealType>(16), static_cast<RealType>(0.25)),
- static_cast<RealType>(0.01)), // Less than cdf == pdf(binomial_distribution<RealType>(16, 0.25), 0)
- static_cast<RealType>(0) // so expect zero as best approximation.
- );
- BOOST_CHECK_EQUAL(
- cdf(
- binomial_distribution<RealType>(static_cast<RealType>(8), static_cast<RealType>(0.25)),
- static_cast<RealType>(8)), static_cast<RealType>(1)
- );
- BOOST_CHECK_EQUAL(
- cdf(
- binomial_distribution<RealType>(static_cast<RealType>(8), static_cast<RealType>(0)),
- static_cast<RealType>(7)), static_cast<RealType>(1)
- );
- BOOST_CHECK_EQUAL(
- cdf(
- binomial_distribution<RealType>(static_cast<RealType>(8), static_cast<RealType>(1)),
- static_cast<RealType>(7)), static_cast<RealType>(0)
- );
- #endif
- {
- // This is a visual sanity check that everything is OK:
- binomial_distribution<RealType> my8dist(8., 0.25); // Note: double values (matching the distribution definition) avoid the need for any casting.
- //cout << "mean(my8dist) = " << boost::math::mean(my8dist) << endl; // mean(my8dist) = 2
- //cout << "my8dist.trials() = " << my8dist.trials() << endl; // my8dist.trials() = 8
- //cout << "my8dist.success_fraction() = " << my8dist.success_fraction() << endl; // my8dist.success_fraction() = 0.25
- BOOST_CHECK_CLOSE(my8dist.trials(), static_cast<RealType>(8), tol2);
- BOOST_CHECK_CLOSE(my8dist.success_fraction(), static_cast<RealType>(0.25), tol2);
- //{
- // int n = static_cast<int>(boost::math::tools::real_cast<double>(my8dist.trials()));
- // RealType sumcdf = 0.;
- // for (int k = 0; k <= n; k++)
- // {
- // cout << k << ' ' << pdf(my8dist, static_cast<RealType>(k));
- // sumcdf += pdf(my8dist, static_cast<RealType>(k));
- // cout << ' ' << sumcdf;
- // cout << ' ' << cdf(my8dist, static_cast<RealType>(k));
- // cout << ' ' << sumcdf - cdf(my8dist, static_cast<RealType>(k)) << endl;
- // } // for k
- // }
- // n = 8, p =0.25
- //k pdf cdf
- //0 0.1001129150390625 0.1001129150390625
- //1 0.26696777343749994 0.36708068847656244
- //2 0.31146240234375017 0.67854309082031261
- //3 0.20764160156249989 0.8861846923828125
- //4 0.086517333984375 0.9727020263671875
- //5 0.023071289062499997 0.9957733154296875
- //6 0.0038452148437500009 0.9996185302734375
- //7 0.00036621093749999984 0.9999847412109375
- //8 1.52587890625e-005 1 1 0
- }
- #define T RealType
- #include "binomial_quantile.ipp"
- for(unsigned i = 0; i < binomial_quantile_data.size(); ++i)
- {
- using namespace boost::math::policies;
- RealType tol = boost::math::tools::epsilon<RealType>() * 500;
- if(!boost::is_floating_point<RealType>::value)
- tol *= 10; // no lanczos approximation implies less accuracy
- RealType x;
- #if !defined(TEST_ROUNDING) || (TEST_ROUNDING == 1)
- //
- // Check full real value first:
- //
- typedef policy<discrete_quantile<boost::math::policies::real> > P1;
- binomial_distribution<RealType, P1> p1(binomial_quantile_data[i][0], binomial_quantile_data[i][1]);
- x = quantile(p1, binomial_quantile_data[i][2]);
- BOOST_CHECK_CLOSE_FRACTION(x, (RealType)binomial_quantile_data[i][3], tol);
- x = quantile(complement(p1, (RealType)binomial_quantile_data[i][2]));
- BOOST_CHECK_CLOSE_FRACTION(x, (RealType)binomial_quantile_data[i][4], tol);
- #endif
- #if !defined(TEST_ROUNDING) || (TEST_ROUNDING == 2)
- //
- // Now with round down to integer:
- //
- typedef policy<discrete_quantile<integer_round_down> > P2;
- binomial_distribution<RealType, P2> p2(binomial_quantile_data[i][0], binomial_quantile_data[i][1]);
- x = quantile(p2, binomial_quantile_data[i][2]);
- BOOST_CHECK_EQUAL(x, (RealType)floor(binomial_quantile_data[i][3]));
- x = quantile(complement(p2, binomial_quantile_data[i][2]));
- BOOST_CHECK_EQUAL(x, (RealType)floor(binomial_quantile_data[i][4]));
- #endif
- #if !defined(TEST_ROUNDING) || (TEST_ROUNDING == 3)
- //
- // Now with round up to integer:
- //
- typedef policy<discrete_quantile<integer_round_up> > P3;
- binomial_distribution<RealType, P3> p3(binomial_quantile_data[i][0], binomial_quantile_data[i][1]);
- x = quantile(p3, binomial_quantile_data[i][2]);
- BOOST_CHECK_EQUAL(x, (RealType)ceil(binomial_quantile_data[i][3]));
- x = quantile(complement(p3, binomial_quantile_data[i][2]));
- BOOST_CHECK_EQUAL(x, (RealType)ceil(binomial_quantile_data[i][4]));
- #endif
- #if !defined(TEST_ROUNDING) || (TEST_ROUNDING == 4)
- //
- // Now with round to integer "outside":
- //
- typedef policy<discrete_quantile<integer_round_outwards> > P4;
- binomial_distribution<RealType, P4> p4(binomial_quantile_data[i][0], binomial_quantile_data[i][1]);
- x = quantile(p4, binomial_quantile_data[i][2]);
- BOOST_CHECK_EQUAL(x, (RealType)(binomial_quantile_data[i][2] < 0.5f ? floor(binomial_quantile_data[i][3]) : ceil(binomial_quantile_data[i][3])));
- x = quantile(complement(p4, binomial_quantile_data[i][2]));
- BOOST_CHECK_EQUAL(x, (RealType)(binomial_quantile_data[i][2] < 0.5f ? ceil(binomial_quantile_data[i][4]) : floor(binomial_quantile_data[i][4])));
- #endif
- #if !defined(TEST_ROUNDING) || (TEST_ROUNDING == 5)
- //
- // Now with round to integer "inside":
- //
- typedef policy<discrete_quantile<integer_round_inwards> > P5;
- binomial_distribution<RealType, P5> p5(binomial_quantile_data[i][0], binomial_quantile_data[i][1]);
- x = quantile(p5, binomial_quantile_data[i][2]);
- BOOST_CHECK_EQUAL(x, (RealType)(binomial_quantile_data[i][2] < 0.5f ? ceil(binomial_quantile_data[i][3]) : floor(binomial_quantile_data[i][3])));
- x = quantile(complement(p5, binomial_quantile_data[i][2]));
- BOOST_CHECK_EQUAL(x, (RealType)(binomial_quantile_data[i][2] < 0.5f ? floor(binomial_quantile_data[i][4]) : ceil(binomial_quantile_data[i][4])));
- #endif
- #if !defined(TEST_ROUNDING) || (TEST_ROUNDING == 6)
- //
- // Now with round to nearest integer:
- //
- typedef policy<discrete_quantile<integer_round_nearest> > P6;
- binomial_distribution<RealType, P6> p6(binomial_quantile_data[i][0], binomial_quantile_data[i][1]);
- x = quantile(p6, binomial_quantile_data[i][2]);
- BOOST_CHECK_EQUAL(x, (RealType)(floor(binomial_quantile_data[i][3] + 0.5f)));
- x = quantile(complement(p6, binomial_quantile_data[i][2]));
- BOOST_CHECK_EQUAL(x, (RealType)(floor(binomial_quantile_data[i][4] + 0.5f)));
- #endif
- }
- check_out_of_range<boost::math::binomial_distribution<RealType> >(1, 1); // (All) valid constructor parameter values.
- } // template <class RealType>void test_spots(RealType)
- BOOST_AUTO_TEST_CASE( test_main )
- {
- BOOST_MATH_CONTROL_FP;
- // Check that can generate binomial distribution using one convenience methods:
- binomial_distribution<> mybn2(1., 0.5); // Using default RealType double.
- // but that
- // boost::math::binomial mybn1(1., 0.5); // Using typedef fails
- // error C2039: 'binomial' : is not a member of 'boost::math'
- // Basic sanity-check spot values.
- // (Parameter value, arbitrarily zero, only communicates the floating point type).
- #ifdef TEST_FLOAT
- test_spots(0.0F); // Test float.
- #endif
- #ifdef TEST_DOUBLE
- test_spots(0.0); // Test double.
- #endif
- #ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
- #ifdef TEST_LDOUBLE
- test_spots(0.0L); // Test long double.
- #endif
- #if !defined(BOOST_MATH_NO_REAL_CONCEPT_TESTS)
- #ifdef TEST_REAL_CONCEPT
- test_spots(boost::math::concepts::real_concept(0.)); // Test real concept.
- #endif
- #endif
- #else
- std::cout << "<note>The long double tests have been disabled on this platform "
- "either because the long double overloads of the usual math functions are "
- "not available at all, or because they are too inaccurate for these tests "
- "to pass.</note>" << std::endl;
- #endif
- } // BOOST_AUTO_TEST_CASE( test_main )
- /*
- Output is:
- Description: Autorun "J:\Cpp\MathToolkit\test\Math_test\Debug\test_binomial.exe"
- Running 1 test case...
- Tolerance for type float is 0.0119209 %
- Tolerance for type double is 2.22045e-011 %
- Tolerance for type long double is 2.22045e-011 %
- Tolerance for type class boost::math::concepts::real_concept is 2.22045e-011 %
- *** No errors detected
- ========== Build: 1 succeeded, 0 failed, 0 up-to-date, 0 skipped ==========
- */
|