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