]>
git.proxmox.com Git - ceph.git/blob - ceph/src/boost/libs/math/test/test_binomial.cpp
3 // Copyright John Maddock 2006.
4 // Copyright Paul A. Bristow 2007.
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)
11 // Basic sanity test for Binomial Cumulative Distribution Function.
13 #define BOOST_MATH_DISCRETE_QUANTILE_POLICY real
15 #if !defined(TEST_FLOAT) && !defined(TEST_DOUBLE) && !defined(TEST_LDOUBLE) && !defined(TEST_REAL_CONCEPT)
19 # define TEST_REAL_CONCEPT
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.
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
;
34 #include <boost/math/distributions/binomial.hpp> // for binomial_distribution
35 using boost::math::binomial_distribution
;
37 #define BOOST_TEST_MAIN
38 #include <boost/test/unit_test.hpp> // for test_main
39 #include <boost/test/floating_point_comparison.hpp> // for BOOST_CHECK_CLOSE
40 #include "table_type.hpp"
42 #include "test_out_of_range.hpp"
48 using std::numeric_limits
;
50 template <class RealType
>
52 RealType N
, // Number of trials
53 RealType k
, // Number of successes
54 RealType p
, // Probability of success
56 RealType Q
, // Complement of CDF
57 RealType tol
) // Test tolerance
59 boost::math::binomial_distribution
<RealType
> bn(N
, p
);
62 if((P
< 0.99) && (Q
< 0.99))
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:
69 cdf(complement(bn
, k
)), Q
, tol
);
73 quantile(bn
, P
), k
, tol
);
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
))
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);
88 quantile(complement(bn
, Q
)), k
, tol
);
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
))
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);
102 // estimate success ratio:
103 // Note lower bound uses a different formual internally
104 // from upper bound, have to adjust things to prevent
107 binomial_distribution
<RealType
>::find_lower_bound_on_p(
111 binomial_distribution
<RealType
>::find_upper_bound_on_p(
117 // Default method (Clopper Pearson)
119 binomial_distribution
<RealType
>::find_lower_bound_on_p(
122 binomial_distribution
<RealType
>::find_upper_bound_on_p(
126 binomial_distribution
<RealType
>::find_lower_bound_on_p(
129 binomial_distribution
<RealType
>::find_upper_bound_on_p(
132 // Bayes Method (Jeffreys Prior)
134 binomial_distribution
<RealType
>::find_lower_bound_on_p(
135 N
, k
, Q
, binomial_distribution
<RealType
>::jeffreys_prior_interval
)
137 binomial_distribution
<RealType
>::find_upper_bound_on_p(
138 N
, k
, Q
, binomial_distribution
<RealType
>::jeffreys_prior_interval
)
141 binomial_distribution
<RealType
>::find_lower_bound_on_p(
142 N
, k
, Q
, binomial_distribution
<RealType
>::jeffreys_prior_interval
)
144 binomial_distribution
<RealType
>::find_upper_bound_on_p(
145 N
, k
, Q
, binomial_distribution
<RealType
>::jeffreys_prior_interval
))
150 // Default method (Clopper Pearson)
152 binomial_distribution
<RealType
>::find_lower_bound_on_p(
155 binomial_distribution
<RealType
>::find_upper_bound_on_p(
159 (binomial_distribution
<RealType
>::find_lower_bound_on_p(
162 binomial_distribution
<RealType
>::find_upper_bound_on_p(
165 // Bayes Method (Jeffreys Prior)
167 binomial_distribution
<RealType
>::find_lower_bound_on_p(
168 N
, k
, P
, binomial_distribution
<RealType
>::jeffreys_prior_interval
)
170 binomial_distribution
<RealType
>::find_upper_bound_on_p(
171 N
, k
, P
, binomial_distribution
<RealType
>::jeffreys_prior_interval
)
174 (binomial_distribution
<RealType
>::find_lower_bound_on_p(
175 N
, k
, P
, binomial_distribution
<RealType
>::jeffreys_prior_interval
)
177 binomial_distribution
<RealType
>::find_upper_bound_on_p(
178 N
, k
, P
, binomial_distribution
<RealType
>::jeffreys_prior_interval
))
183 // estimate sample size:
186 binomial_distribution
<RealType
>::find_minimum_number_of_trials(
190 binomial_distribution
<RealType
>::find_maximum_number_of_trials(
195 // Double check consistency of CDF and PDF by computing
198 for(unsigned i
= 0; i
<= k
; ++i
)
199 sum
+= pdf(bn
, RealType(i
));
202 // And complement as well:
204 for(RealType i
= N
; i
> k
; i
-= 1)
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
);
219 template <class RealType
> // Any floating-point type RealType.
220 void test_spots(RealType T
)
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
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
233 cout
<< "Tolerance for type " << typeid(T
).name() << " is " << tolerance
<< " %" << endl
;
236 // Sources of spot test values:
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.
242 // P = pbinom(30, 500, 0.05) = 0.869147702104609
244 using boost::math::binomial_distribution
;
245 using ::boost::math::cdf
;
246 using ::boost::math::pdf
;
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.
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
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
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!
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
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
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
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
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
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
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
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
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
347 // Tests on PDF follow:
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
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
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
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
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
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
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
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
416 // Some exact (probably) values.
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
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
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
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
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
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
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.
459 , static_cast<RealType
>(8 * 0.25), tol2
);
463 , static_cast<RealType
>(8 * 0.25 * 0.75), tol2
);
466 standard_deviation(dist
)
467 , static_cast<RealType
>(sqrt(8 * 0.25L * 0.75L)), tol2
);
471 , pdf(dist
, x
) / cdf(complement(dist
, x
)), tol2
);
472 // cumulative hazard:
475 , -log(cdf(complement(dist
, x
))), tol2
);
476 // coefficient_of_variation:
478 coefficient_of_variation(dist
)
479 , standard_deviation(dist
) / mean(dist
), tol2
);
483 , static_cast<RealType
>(std::floor(9 * 0.25)), tol2
);
487 , static_cast<RealType
>(0.40824829046386301636621401245098L), (std::max
)(tol2
, static_cast<RealType
>(5e-29))); // test data has 32 digits only.
491 , static_cast<RealType
>(2.916666666666666666666666666666666666L), tol2
);
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
));
499 // special cases for PDF:
502 binomial_distribution
<RealType
>(static_cast<RealType
>(8), static_cast<RealType
>(0)),
503 static_cast<RealType
>(0)), static_cast<RealType
>(1)
507 binomial_distribution
<RealType
>(static_cast<RealType
>(8), static_cast<RealType
>(0)),
508 static_cast<RealType
>(0.0001)), static_cast<RealType
>(0)
512 binomial_distribution
<RealType
>(static_cast<RealType
>(8), static_cast<RealType
>(1)),
513 static_cast<RealType
>(0.001)), static_cast<RealType
>(0)
517 binomial_distribution
<RealType
>(static_cast<RealType
>(8), static_cast<RealType
>(1)),
518 static_cast<RealType
>(8)), static_cast<RealType
>(1)
522 binomial_distribution
<RealType
>(static_cast<RealType
>(0), static_cast<RealType
>(0.25)),
523 static_cast<RealType
>(0)), static_cast<RealType
>(1)
525 BOOST_MATH_CHECK_THROW(
527 binomial_distribution
<RealType
>(static_cast<RealType
>(-1), static_cast<RealType
>(0.25)),
528 static_cast<RealType
>(0)), std::domain_error
530 BOOST_MATH_CHECK_THROW(
532 binomial_distribution
<RealType
>(static_cast<RealType
>(8), static_cast<RealType
>(-0.25)),
533 static_cast<RealType
>(0)), std::domain_error
535 BOOST_MATH_CHECK_THROW(
537 binomial_distribution
<RealType
>(static_cast<RealType
>(8), static_cast<RealType
>(1.25)),
538 static_cast<RealType
>(0)), std::domain_error
540 BOOST_MATH_CHECK_THROW(
542 binomial_distribution
<RealType
>(static_cast<RealType
>(8), static_cast<RealType
>(0.25)),
543 static_cast<RealType
>(-1)), std::domain_error
545 BOOST_MATH_CHECK_THROW(
547 binomial_distribution
<RealType
>(static_cast<RealType
>(8), static_cast<RealType
>(0.25)),
548 static_cast<RealType
>(9)), std::domain_error
550 BOOST_MATH_CHECK_THROW(
552 binomial_distribution
<RealType
>(static_cast<RealType
>(8), static_cast<RealType
>(0.25)),
553 static_cast<RealType
>(-1)), std::domain_error
555 BOOST_MATH_CHECK_THROW(
557 binomial_distribution
<RealType
>(static_cast<RealType
>(8), static_cast<RealType
>(0.25)),
558 static_cast<RealType
>(9)), std::domain_error
560 BOOST_MATH_CHECK_THROW(
562 binomial_distribution
<RealType
>(static_cast<RealType
>(8), static_cast<RealType
>(-0.25)),
563 static_cast<RealType
>(0)), std::domain_error
565 BOOST_MATH_CHECK_THROW(
567 binomial_distribution
<RealType
>(static_cast<RealType
>(8), static_cast<RealType
>(1.25)),
568 static_cast<RealType
>(0)), std::domain_error
570 BOOST_MATH_CHECK_THROW(
572 binomial_distribution
<RealType
>(static_cast<RealType
>(8), static_cast<RealType
>(-0.25)),
573 static_cast<RealType
>(0)), std::domain_error
575 BOOST_MATH_CHECK_THROW(
577 binomial_distribution
<RealType
>(static_cast<RealType
>(8), static_cast<RealType
>(1.25)),
578 static_cast<RealType
>(0)), std::domain_error
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.
590 binomial_distribution
<RealType
>(static_cast<RealType
>(8), static_cast<RealType
>(0.25)),
591 static_cast<RealType
>(8)), static_cast<RealType
>(1)
595 binomial_distribution
<RealType
>(static_cast<RealType
>(8), static_cast<RealType
>(0)),
596 static_cast<RealType
>(7)), static_cast<RealType
>(1)
600 binomial_distribution
<RealType
>(static_cast<RealType
>(8), static_cast<RealType
>(1)),
601 static_cast<RealType
>(7)), static_cast<RealType
>(0)
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
);
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++)
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;
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
640 #include "binomial_quantile.ipp"
642 for(unsigned i
= 0; i
< binomial_quantile_data
.size(); ++i
)
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
649 #if !defined(TEST_ROUNDING) || (TEST_ROUNDING == 1)
651 // Check full real value first:
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
);
660 #if !defined(TEST_ROUNDING) || (TEST_ROUNDING == 2)
662 // Now with round down to integer:
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]));
671 #if !defined(TEST_ROUNDING) || (TEST_ROUNDING == 3)
673 // Now with round up to integer:
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]));
682 #if !defined(TEST_ROUNDING) || (TEST_ROUNDING == 4)
684 // Now with round to integer "outside":
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])));
693 #if !defined(TEST_ROUNDING) || (TEST_ROUNDING == 5)
695 // Now with round to integer "inside":
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])));
704 #if !defined(TEST_ROUNDING) || (TEST_ROUNDING == 6)
706 // Now with round to nearest integer:
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
)));
717 check_out_of_range
<boost::math::binomial_distribution
<RealType
> >(1, 1); // (All) valid constructor parameter values.
720 } // template <class RealType>void test_spots(RealType)
722 BOOST_AUTO_TEST_CASE( test_main
)
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.
728 // boost::math::binomial mybn1(1., 0.5); // Using typedef fails
729 // error C2039: 'binomial' : is not a member of 'boost::math'
731 // Basic sanity-check spot values.
733 // (Parameter value, arbitrarily zero, only communicates the floating point type).
735 test_spots(0.0F
); // Test float.
738 test_spots(0.0); // Test double.
740 #ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
742 test_spots(0.0L); // Test long double.
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.
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
;
756 } // BOOST_AUTO_TEST_CASE( test_main )
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 %
769 *** No errors detected
771 ========== Build: 1 succeeded, 0 failed, 0 up-to-date, 0 skipped ==========