1 // test_negative_binomial.cpp
2
3 // Copyright Paul A. Bristow 2007.
4 // Copyright John Maddock 2006.
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 // Tests for Negative Binomial Distribution.
12
13 // Note that these defines must be placed BEFORE #includes.
14 #define BOOST_MATH_OVERFLOW_ERROR_POLICY ignore_error
15 // because several tests overflow & underflow by design.
16 #define BOOST_MATH_DISCRETE_QUANTILE_POLICY real
17
18 #ifdef _MSC_VER
19 # pragma warning(disable: 4127) // conditional expression is constant.
20 #endif
21
22 #if !defined(TEST_FLOAT) && !defined(TEST_DOUBLE) && !defined(TEST_LDOUBLE) && !defined(TEST_REAL_CONCEPT)
23 # define TEST_FLOAT
24 # define TEST_DOUBLE
25 # define TEST_LDOUBLE
26 # define TEST_REAL_CONCEPT
27 #endif
28
29 #include <boost/math/concepts/real_concept.hpp> // for real_concept
30 using ::boost::math::concepts::real_concept;
31
32 #include <boost/math/distributions/negative_binomial.hpp> // for negative_binomial_distribution
33 using boost::math::negative_binomial_distribution;
34
35 #include <boost/math/special_functions/gamma.hpp>
36 using boost::math::lgamma; // log gamma
37
38 #define BOOST_TEST_MAIN
39 #include <boost/test/unit_test.hpp> // for test_main
40 #include <boost/test/floating_point_comparison.hpp> // for BOOST_CHECK_CLOSE
41 #include "table_type.hpp"
42 #include "test_out_of_range.hpp"
43
44 #include <iostream>
45 using std::cout;
46 using std::endl;
47 using std::setprecision;
48 using std::showpoint;
49 #include <limits>
50 using std::numeric_limits;
51
52 template <class RealType>
test_spot(RealType N,RealType k,RealType p,RealType P,RealType Q,RealType tol)53 void test_spot( // Test a single spot value against 'known good' values.
54 RealType N, // Number of successes.
55 RealType k, // Number of failures.
56 RealType p, // Probability of success_fraction.
57 RealType P, // CDF probability.
58 RealType Q, // Complement of CDF.
59 RealType tol) // Test tolerance.
60 {
61 boost::math::negative_binomial_distribution<RealType> bn(N, p);
62 BOOST_CHECK_EQUAL(N, bn.successes());
63 BOOST_CHECK_EQUAL(p, bn.success_fraction());
64 BOOST_CHECK_CLOSE(
65 cdf(bn, k), P, tol);
66
67 if((P < 0.99) && (Q < 0.99))
68 {
69 // We can only check this if P is not too close to 1,
70 // so that we can guarantee that Q is free of error:
71 //
72 BOOST_CHECK_CLOSE(
73 cdf(complement(bn, k)), Q, tol);
74 if(k != 0)
75 {
76 BOOST_CHECK_CLOSE(
77 quantile(bn, P), k, tol);
78 }
79 else
80 {
81 // Just check quantile is very small:
82 if((std::numeric_limits<RealType>::max_exponent <= std::numeric_limits<double>::max_exponent)
83 && (boost::is_floating_point<RealType>::value))
84 {
85 // Limit where this is checked: if exponent range is very large we may
86 // run out of iterations in our root finding algorithm.
87 BOOST_CHECK(quantile(bn, P) < boost::math::tools::epsilon<RealType>() * 10);
88 }
89 }
90 if(k != 0)
91 {
92 BOOST_CHECK_CLOSE(
93 quantile(complement(bn, Q)), k, tol);
94 }
95 else
96 {
97 // Just check quantile is very small:
98 if((std::numeric_limits<RealType>::max_exponent <= std::numeric_limits<double>::max_exponent)
99 && (boost::is_floating_point<RealType>::value))
100 {
101 // Limit where this is checked: if exponent range is very large we may
102 // run out of iterations in our root finding algorithm.
103 BOOST_CHECK(quantile(complement(bn, Q)) < boost::math::tools::epsilon<RealType>() * 10);
104 }
105 }
106 // estimate success ratio:
107 BOOST_CHECK_CLOSE(
108 negative_binomial_distribution<RealType>::find_lower_bound_on_p(
109 N+k, N, P),
110 p, tol);
111 // Note we bump up the sample size here, purely for the sake of the test,
112 // internally the function has to adjust the sample size so that we get
113 // the right upper bound, our test undoes this, so we can verify the result.
114 BOOST_CHECK_CLOSE(
115 negative_binomial_distribution<RealType>::find_upper_bound_on_p(
116 N+k+1, N, Q),
117 p, tol);
118
119 if(Q < P)
120 {
121 //
122 // We check two things here, that the upper and lower bounds
123 // are the right way around, and that they do actually bracket
124 // the naive estimate of p = successes / (sample size)
125 //
126 BOOST_CHECK(
127 negative_binomial_distribution<RealType>::find_lower_bound_on_p(
128 N+k, N, Q)
129 <=
130 negative_binomial_distribution<RealType>::find_upper_bound_on_p(
131 N+k, N, Q)
132 );
133 BOOST_CHECK(
134 negative_binomial_distribution<RealType>::find_lower_bound_on_p(
135 N+k, N, Q)
136 <=
137 N / (N+k)
138 );
139 BOOST_CHECK(
140 N / (N+k)
141 <=
142 negative_binomial_distribution<RealType>::find_upper_bound_on_p(
143 N+k, N, Q)
144 );
145 }
146 else
147 {
148 // As above but when P is small.
149 BOOST_CHECK(
150 negative_binomial_distribution<RealType>::find_lower_bound_on_p(
151 N+k, N, P)
152 <=
153 negative_binomial_distribution<RealType>::find_upper_bound_on_p(
154 N+k, N, P)
155 );
156 BOOST_CHECK(
157 negative_binomial_distribution<RealType>::find_lower_bound_on_p(
158 N+k, N, P)
159 <=
160 N / (N+k)
161 );
162 BOOST_CHECK(
163 N / (N+k)
164 <=
165 negative_binomial_distribution<RealType>::find_upper_bound_on_p(
166 N+k, N, P)
167 );
168 }
169
170 // Estimate sample size:
171 BOOST_CHECK_CLOSE(
172 negative_binomial_distribution<RealType>::find_minimum_number_of_trials(
173 k, p, P),
174 N+k, tol);
175 BOOST_CHECK_CLOSE(
176 negative_binomial_distribution<RealType>::find_maximum_number_of_trials(
177 k, p, Q),
178 N+k, tol);
179
180 // Double check consistency of CDF and PDF by computing the finite sum:
181 RealType sum = 0;
182 for(unsigned i = 0; i <= k; ++i)
183 {
184 sum += pdf(bn, RealType(i));
185 }
186 BOOST_CHECK_CLOSE(sum, P, tol);
187
188 // Complement is not possible since sum is to infinity.
189 } //
190 } // test_spot
191
192 template <class RealType> // Any floating-point type RealType.
test_spots(RealType)193 void test_spots(RealType)
194 {
195 // Basic sanity checks, test data is to double precision only
196 // so set tolerance to 1000 eps expressed as a percent, or
197 // 1000 eps of type double expressed as a percent, whichever
198 // is the larger.
199
200 RealType tolerance = (std::max)
201 (boost::math::tools::epsilon<RealType>(),
202 static_cast<RealType>(std::numeric_limits<double>::epsilon()));
203 tolerance *= 100 * 100000.0f;
204
205 cout << "Tolerance = " << tolerance << "%." << endl;
206
207 RealType tol1eps = boost::math::tools::epsilon<RealType>() * 2; // Very tight, suit exact values.
208 //RealType tol2eps = boost::math::tools::epsilon<RealType>() * 2; // Tight, suit exact values.
209 RealType tol5eps = boost::math::tools::epsilon<RealType>() * 5; // Wider 5 epsilon.
210 cout << "Tolerance 5 eps = " << tol5eps << "%." << endl;
211
212 // Sources of spot test values:
213
214 // MathCAD defines pbinom(k, r, p) (at about 64-bit double precision, about 16 decimal digits)
215 // returns pr(X , k) when random variable X has the binomial distribution with parameters r and p.
216 // 0 <= k
217 // r > 0
218 // 0 <= p <= 1
219 // P = pbinom(30, 500, 0.05) = 0.869147702104609
220
221 // And functions.wolfram.com
222
223 using boost::math::negative_binomial_distribution;
224 using ::boost::math::negative_binomial;
225 using ::boost::math::cdf;
226 using ::boost::math::pdf;
227
228 // Test negative binomial using cdf spot values from MathCAD cdf = pnbinom(k, r, p).
229 // These test quantiles and complements as well.
230
231 test_spot( // pnbinom(1,2,0.5) = 0.5
232 static_cast<RealType>(2), // successes r
233 static_cast<RealType>(1), // Number of failures, k
234 static_cast<RealType>(0.5), // Probability of success as fraction, p
235 static_cast<RealType>(0.5), // Probability of result (CDF), P
236 static_cast<RealType>(0.5), // complement CCDF Q = 1 - P
237 tolerance);
238
239 test_spot( // pbinom(0, 2, 0.25)
240 static_cast<RealType>(2), // successes r
241 static_cast<RealType>(0), // Number of failures, k
242 static_cast<RealType>(0.25),
243 static_cast<RealType>(0.0625), // Probability of result (CDF), P
244 static_cast<RealType>(0.9375), // Q = 1 - P
245 tolerance);
246
247 test_spot( // pbinom(48,8,0.25)
248 static_cast<RealType>(8), // successes r
249 static_cast<RealType>(48), // Number of failures, k
250 static_cast<RealType>(0.25), // Probability of success, p
251 static_cast<RealType>(9.826582228110670E-1), // Probability of result (CDF), P
252 static_cast<RealType>(1 - 9.826582228110670E-1), // Q = 1 - P
253 tolerance);
254
255 test_spot( // pbinom(2,5,0.4)
256 static_cast<RealType>(5), // successes r
257 static_cast<RealType>(2), // Number of failures, k
258 static_cast<RealType>(0.4), // Probability of success, p
259 static_cast<RealType>(9.625600000000020E-2), // Probability of result (CDF), P
260 static_cast<RealType>(1 - 9.625600000000020E-2), // Q = 1 - P
261 tolerance);
262
263 test_spot( // pbinom(10,100,0.9)
264 static_cast<RealType>(100), // successes r
265 static_cast<RealType>(10), // Number of failures, k
266 static_cast<RealType>(0.9), // Probability of success, p
267 static_cast<RealType>(4.535522887695670E-1), // Probability of result (CDF), P
268 static_cast<RealType>(1 - 4.535522887695670E-1), // Q = 1 - P
269 tolerance);
270
271 test_spot( // pbinom(1,100,0.991)
272 static_cast<RealType>(100), // successes r
273 static_cast<RealType>(1), // Number of failures, k
274 static_cast<RealType>(0.991), // Probability of success, p
275 static_cast<RealType>(7.693413044217000E-1), // Probability of result (CDF), P
276 static_cast<RealType>(1 - 7.693413044217000E-1), // Q = 1 - P
277 tolerance);
278
279 test_spot( // pbinom(10,100,0.991)
280 static_cast<RealType>(100), // successes r
281 static_cast<RealType>(10), // Number of failures, k
282 static_cast<RealType>(0.991), // Probability of success, p
283 static_cast<RealType>(9.999999940939000E-1), // Probability of result (CDF), P
284 static_cast<RealType>(1 - 9.999999940939000E-1), // Q = 1 - P
285 tolerance);
286
287 if(std::numeric_limits<RealType>::is_specialized)
288 { // An extreme value test that takes 3 minutes using the real concept type
289 // for which numeric_limits<RealType>::is_specialized == false, deliberately
290 // and for which there is no Lanczos approximation defined (also deliberately)
291 // giving a very slow computation, but with acceptable accuracy.
292 // A possible enhancement might be to use a normal approximation for
293 // extreme values, but this is not implemented.
294 test_spot( // pbinom(100000,100,0.001)
295 static_cast<RealType>(100), // successes r
296 static_cast<RealType>(100000), // Number of failures, k
297 static_cast<RealType>(0.001), // Probability of success, p
298 static_cast<RealType>(5.173047534260320E-1), // Probability of result (CDF), P
299 static_cast<RealType>(1 - 5.173047534260320E-1), // Q = 1 - P
300 tolerance*1000); // *1000 is OK 0.51730475350664229 versus
301
302 // functions.wolfram.com
303 // for I[0.001](100, 100000+1) gives:
304 // Wolfram 0.517304753506834882009032744488738352004003696396461766326713
305 // JM nonLanczos 0.51730475350664229 differs at the 13th decimal digit.
306 // MathCAD 0.51730475342603199 differs at 10th decimal digit.
307
308 // Error tests:
309 check_out_of_range<negative_binomial_distribution<RealType> >(20, 0.5);
310 BOOST_CHECK_THROW(negative_binomial_distribution<RealType>(0, 0.5), std::domain_error);
311 BOOST_CHECK_THROW(negative_binomial_distribution<RealType>(-2, 0.5), std::domain_error);
312 BOOST_CHECK_THROW(negative_binomial_distribution<RealType>(20, -0.5), std::domain_error);
313 BOOST_CHECK_THROW(negative_binomial_distribution<RealType>(20, 1.5), std::domain_error);
314 }
315 // End of single spot tests using RealType
316
317
318 // Tests on PDF:
319 BOOST_CHECK_CLOSE(
320 pdf(negative_binomial_distribution<RealType>(static_cast<RealType>(2), static_cast<RealType>(0.5)),
321 static_cast<RealType>(0) ), // k = 0.
322 static_cast<RealType>(0.25), // 0
323 tolerance);
324
325 BOOST_CHECK_CLOSE(
326 pdf(negative_binomial_distribution<RealType>(static_cast<RealType>(4), static_cast<RealType>(0.5)),
327 static_cast<RealType>(0)), // k = 0.
328 static_cast<RealType>(0.0625), // exact 1/16
329 tolerance);
330
331 BOOST_CHECK_CLOSE(
332 pdf(negative_binomial_distribution<RealType>(static_cast<RealType>(20), static_cast<RealType>(0.25)),
333 static_cast<RealType>(0)), // k = 0
334 static_cast<RealType>(9.094947017729270E-13), // pbinom(0,20,0.25) = 9.094947017729270E-13
335 tolerance);
336
337 BOOST_CHECK_CLOSE(
338 pdf(negative_binomial_distribution<RealType>(static_cast<RealType>(20), static_cast<RealType>(0.2)),
339 static_cast<RealType>(0)), // k = 0
340 static_cast<RealType>(1.0485760000000003e-014), // MathCAD 1.048576000000000E-14
341 tolerance);
342
343 BOOST_CHECK_CLOSE(
344 pdf(negative_binomial_distribution<RealType>(static_cast<RealType>(10), static_cast<RealType>(0.1)),
345 static_cast<RealType>(0)), // k = 0.
346 static_cast<RealType>(1e-10), // MathCAD says zero, but suffers cancellation error?
347 tolerance);
348
349 BOOST_CHECK_CLOSE(
350 pdf(negative_binomial_distribution<RealType>(static_cast<RealType>(20), static_cast<RealType>(0.1)),
351 static_cast<RealType>(0)), // k = 0.
352 static_cast<RealType>(1e-20), // MathCAD says zero, but suffers cancellation error?
353 tolerance);
354
355
356 BOOST_CHECK_CLOSE( // .
357 pdf(negative_binomial_distribution<RealType>(static_cast<RealType>(20), static_cast<RealType>(0.9)),
358 static_cast<RealType>(0)), // k.
359 static_cast<RealType>(1.215766545905690E-1), // k=20 p = 0.9
360 tolerance);
361
362 // Tests on cdf:
363 // MathCAD pbinom k, r, p) == failures, successes, probability.
364
365 BOOST_CHECK_CLOSE(cdf(
366 negative_binomial_distribution<RealType>(static_cast<RealType>(2), static_cast<RealType>(0.5)), // successes = 2,prob 0.25
367 static_cast<RealType>(0) ), // k = 0
368 static_cast<RealType>(0.25), // probability 1/4
369 tolerance);
370
371 BOOST_CHECK_CLOSE(cdf(complement(
372 negative_binomial_distribution<RealType>(static_cast<RealType>(2), static_cast<RealType>(0.5)), // successes = 2,prob 0.25
373 static_cast<RealType>(0) )), // k = 0
374 static_cast<RealType>(0.75), // probability 3/4
375 tolerance);
376 BOOST_CHECK_CLOSE( // k = 1.
377 cdf(negative_binomial_distribution<RealType>(static_cast<RealType>(20), static_cast<RealType>(0.25)),
378 static_cast<RealType>(1)), // k =1.
379 static_cast<RealType>(1.455191522836700E-11),
380 tolerance);
381
382 BOOST_CHECK_SMALL( // Check within an epsilon with CHECK_SMALL
383 cdf(negative_binomial_distribution<RealType>(static_cast<RealType>(20), static_cast<RealType>(0.25)),
384 static_cast<RealType>(1)) -
385 static_cast<RealType>(1.455191522836700E-11),
386 tolerance );
387
388 // Some exact (probably - judging by trailing zeros) values.
389 BOOST_CHECK_CLOSE(
390 cdf(negative_binomial_distribution<RealType>(static_cast<RealType>(8), static_cast<RealType>(0.25)),
391 static_cast<RealType>(0)), // k.
392 static_cast<RealType>(1.525878906250000E-5),
393 tolerance);
394
395 BOOST_CHECK_CLOSE(
396 cdf(negative_binomial_distribution<RealType>(static_cast<RealType>(8), static_cast<RealType>(0.25)),
397 static_cast<RealType>(0)), // k.
398 static_cast<RealType>(1.525878906250000E-5),
399 tolerance);
400
401 BOOST_CHECK_SMALL(
402 cdf(negative_binomial_distribution<RealType>(static_cast<RealType>(8), static_cast<RealType>(0.25)),
403 static_cast<RealType>(0)) -
404 static_cast<RealType>(1.525878906250000E-5),
405 tolerance );
406
407 BOOST_CHECK_CLOSE( // k = 1.
408 cdf(negative_binomial_distribution<RealType>(static_cast<RealType>(8), static_cast<RealType>(0.25)),
409 static_cast<RealType>(1)), // k.
410 static_cast<RealType>(1.068115234375010E-4),
411 tolerance);
412
413 BOOST_CHECK_CLOSE( // k = 2.
414 cdf(negative_binomial_distribution<RealType>(static_cast<RealType>(8), static_cast<RealType>(0.25)),
415 static_cast<RealType>(2)), // k.
416 static_cast<RealType>(4.158020019531300E-4),
417 tolerance);
418
419 BOOST_CHECK_CLOSE( // k = 3.
420 cdf(negative_binomial_distribution<RealType>(static_cast<RealType>(8), static_cast<RealType>(0.25)),
421 static_cast<RealType>(3)), // k.bristow
422 static_cast<RealType>(1.188278198242200E-3),
423 tolerance);
424
425 BOOST_CHECK_CLOSE( // k = 4.
426 cdf(negative_binomial_distribution<RealType>(static_cast<RealType>(8), static_cast<RealType>(0.25)),
427 static_cast<RealType>(4)), // k.
428 static_cast<RealType>(2.781510353088410E-3),
429 tolerance);
430
431 BOOST_CHECK_CLOSE( // k = 5.
432 cdf(negative_binomial_distribution<RealType>(static_cast<RealType>(8), static_cast<RealType>(0.25)),
433 static_cast<RealType>(5)), // k.
434 static_cast<RealType>(5.649328231811500E-3),
435 tolerance);
436
437 BOOST_CHECK_CLOSE( // k = 6.
438 cdf(negative_binomial_distribution<RealType>(static_cast<RealType>(8), static_cast<RealType>(0.25)),
439 static_cast<RealType>(6)), // k.
440 static_cast<RealType>(1.030953228473680E-2),
441 tolerance);
442
443 BOOST_CHECK_CLOSE( // k = 7.
444 cdf(negative_binomial_distribution<RealType>(static_cast<RealType>(8), static_cast<RealType>(0.25)),
445 static_cast<RealType>(7)), // k.
446 static_cast<RealType>(1.729983836412430E-2),
447 tolerance);
448
449 BOOST_CHECK_CLOSE( // k = 8.
450 cdf(negative_binomial_distribution<RealType>(static_cast<RealType>(8), static_cast<RealType>(0.25)),
451 static_cast<RealType>(8)), // k = n.
452 static_cast<RealType>(2.712995628826370E-2),
453 tolerance);
454
455 BOOST_CHECK_CLOSE( //
456 cdf(negative_binomial_distribution<RealType>(static_cast<RealType>(8), static_cast<RealType>(0.25)),
457 static_cast<RealType>(48)), // k
458 static_cast<RealType>(9.826582228110670E-1),
459 tolerance);
460
461 BOOST_CHECK_CLOSE( //
462 cdf(negative_binomial_distribution<RealType>(static_cast<RealType>(8), static_cast<RealType>(0.25)),
463 static_cast<RealType>(64)), // k
464 static_cast<RealType>(9.990295004935590E-1),
465 tolerance);
466
467 BOOST_CHECK_CLOSE( //
468 cdf(negative_binomial_distribution<RealType>(static_cast<RealType>(5), static_cast<RealType>(0.4)),
469 static_cast<RealType>(26)), // k
470 static_cast<RealType>(9.989686246611190E-1),
471 tolerance);
472
473 BOOST_CHECK_CLOSE( //
474 cdf(negative_binomial_distribution<RealType>(static_cast<RealType>(5), static_cast<RealType>(0.4)),
475 static_cast<RealType>(2)), // k failures
476 static_cast<RealType>(9.625600000000020E-2),
477 tolerance);
478
479 BOOST_CHECK_CLOSE( //
480 cdf(negative_binomial_distribution<RealType>(static_cast<RealType>(50), static_cast<RealType>(0.9)),
481 static_cast<RealType>(20)), // k
482 static_cast<RealType>(9.999970854144170E-1),
483 tolerance);
484
485 BOOST_CHECK_CLOSE( //
486 cdf(negative_binomial_distribution<RealType>(static_cast<RealType>(500), static_cast<RealType>(0.7)),
487 static_cast<RealType>(200)), // k
488 static_cast<RealType>(2.172846379930550E-1),
489 tolerance* 2);
490
491 BOOST_CHECK_CLOSE( //
492 cdf(negative_binomial_distribution<RealType>(static_cast<RealType>(50), static_cast<RealType>(0.7)),
493 static_cast<RealType>(20)), // k
494 static_cast<RealType>(4.550203671301790E-1),
495 tolerance);
496
497 // Tests of other functions, mean and other moments ...
498
499 negative_binomial_distribution<RealType> dist(static_cast<RealType>(8), static_cast<RealType>(0.25));
500 using namespace std; // ADL of std names.
501 // mean:
502 BOOST_CHECK_CLOSE(
503 mean(dist), static_cast<RealType>(8 * (1 - 0.25) /0.25), tol5eps);
504 BOOST_CHECK_CLOSE(
505 mode(dist), static_cast<RealType>(21), tol1eps);
506 // variance:
507 BOOST_CHECK_CLOSE(
508 variance(dist), static_cast<RealType>(8 * (1 - 0.25) / (0.25 * 0.25)), tol5eps);
509 // std deviation:
510 BOOST_CHECK_CLOSE(
511 standard_deviation(dist), // 9.79795897113271239270
512 static_cast<RealType>(9.797958971132712392789136298823565567864L), // using functions.wolfram.com
513 // 9.79795897113271152534 == sqrt(8 * (1 - 0.25) / (0.25 * 0.25)))
514 tol5eps * 100);
515 BOOST_CHECK_CLOSE(
516 skewness(dist), //
517 static_cast<RealType>(0.71443450831176036),
518 // using http://mathworld.wolfram.com/skewness.html
519 tolerance);
520 BOOST_CHECK_CLOSE(
521 kurtosis_excess(dist), //
522 static_cast<RealType>(0.7604166666666666666666666666666666666666L), // using Wikipedia Kurtosis(excess) formula
523 tol5eps * 100);
524 BOOST_CHECK_CLOSE(
525 kurtosis(dist), // true
526 static_cast<RealType>(3.76041666666666666666666666666666666666666L), //
527 tol5eps * 100);
528 // hazard:
529 RealType x = static_cast<RealType>(0.125);
530 BOOST_CHECK_CLOSE(
531 hazard(dist, x)
532 , pdf(dist, x) / cdf(complement(dist, x)), tol5eps);
533 // cumulative hazard:
534 BOOST_CHECK_CLOSE(
535 chf(dist, x), -log(cdf(complement(dist, x))), tol5eps);
536 // coefficient_of_variation:
537 BOOST_CHECK_CLOSE(
538 coefficient_of_variation(dist)
539 , standard_deviation(dist) / mean(dist), tol5eps);
540
541 // Special cases for PDF:
542 BOOST_CHECK_EQUAL(
543 pdf(
544 negative_binomial_distribution<RealType>(static_cast<RealType>(8), static_cast<RealType>(0)), //
545 static_cast<RealType>(0)),
546 static_cast<RealType>(0) );
547
548 BOOST_CHECK_EQUAL(
549 pdf(
550 negative_binomial_distribution<RealType>(static_cast<RealType>(8), static_cast<RealType>(0)),
551 static_cast<RealType>(0.0001)),
552 static_cast<RealType>(0) );
553
554 BOOST_CHECK_EQUAL(
555 pdf(
556 negative_binomial_distribution<RealType>(static_cast<RealType>(8), static_cast<RealType>(1)),
557 static_cast<RealType>(0.001)),
558 static_cast<RealType>(0) );
559
560 BOOST_CHECK_EQUAL(
561 pdf(
562 negative_binomial_distribution<RealType>(static_cast<RealType>(8), static_cast<RealType>(1)),
563 static_cast<RealType>(8)),
564 static_cast<RealType>(0) );
565
566 BOOST_CHECK_SMALL(
567 pdf(
568 negative_binomial_distribution<RealType>(static_cast<RealType>(2), static_cast<RealType>(0.25)),
569 static_cast<RealType>(0))-
570 static_cast<RealType>(0.0625),
571 2 * boost::math::tools::epsilon<RealType>() ); // Expect exact, but not quite.
572 // numeric_limits<RealType>::epsilon()); // Not suitable for real concept!
573
574 // Quantile boundary cases checks:
575 BOOST_CHECK_EQUAL(
576 quantile( // zero P < cdf(0) so should be exactly zero.
577 negative_binomial_distribution<RealType>(static_cast<RealType>(8), static_cast<RealType>(0.25)),
578 static_cast<RealType>(0)),
579 static_cast<RealType>(0));
580
581 BOOST_CHECK_EQUAL(
582 quantile( // min P < cdf(0) so should be exactly zero.
583 negative_binomial_distribution<RealType>(static_cast<RealType>(8), static_cast<RealType>(0.25)),
584 static_cast<RealType>(boost::math::tools::min_value<RealType>())),
585 static_cast<RealType>(0));
586
587 BOOST_CHECK_CLOSE_FRACTION(
588 quantile( // Small P < cdf(0) so should be near zero.
589 negative_binomial_distribution<RealType>(static_cast<RealType>(8), static_cast<RealType>(0.25)),
590 static_cast<RealType>(boost::math::tools::epsilon<RealType>())), //
591 static_cast<RealType>(0),
592 tol5eps);
593
594 BOOST_CHECK_CLOSE(
595 quantile( // Small P < cdf(0) so should be exactly zero.
596 negative_binomial_distribution<RealType>(static_cast<RealType>(8), static_cast<RealType>(0.25)),
597 static_cast<RealType>(0.0001)),
598 static_cast<RealType>(0.95854156929288470),
599 tolerance);
600
601 //BOOST_CHECK( // Fails with overflow for real_concept
602 //quantile( // Small P near 1 so k failures should be big.
603 //negative_binomial_distribution<RealType>(static_cast<RealType>(8), static_cast<RealType>(0.25)),
604 //static_cast<RealType>(1 - boost::math::tools::epsilon<RealType>())) <=
605 //static_cast<RealType>(189.56999032670058) // 106.462769 for float
606 //);
607
608 if(std::numeric_limits<RealType>::has_infinity)
609 { // BOOST_CHECK tests for infinity using std::numeric_limits<>::infinity()
610 // Note that infinity is not implemented for real_concept, so these tests
611 // are only done for types, like built-in float, double.. that have infinity.
612 // Note that these assume that BOOST_MATH_OVERFLOW_ERROR_POLICY is NOT throw_on_error.
613 // #define BOOST_MATH_THROW_ON_OVERFLOW_POLICY == throw_on_error would throw here.
614 // #define BOOST_MAT_DOMAIN_ERROR_POLICY IS defined throw_on_error,
615 // so the throw path of error handling is tested below with BOOST_CHECK_THROW tests.
616
617 BOOST_CHECK(
618 quantile( // At P == 1 so k failures should be infinite.
619 negative_binomial_distribution<RealType>(static_cast<RealType>(8), static_cast<RealType>(0.25)),
620 static_cast<RealType>(1)) ==
621 //static_cast<RealType>(boost::math::tools::infinity<RealType>())
622 static_cast<RealType>(std::numeric_limits<RealType>::infinity()) );
623
624 BOOST_CHECK_EQUAL(
625 quantile( // At 1 == P so should be infinite.
626 negative_binomial_distribution<RealType>(static_cast<RealType>(8), static_cast<RealType>(0.25)),
627 static_cast<RealType>(1)), //
628 std::numeric_limits<RealType>::infinity() );
629
630 BOOST_CHECK_EQUAL(
631 quantile(complement( // Q zero 1 so P == 1 < cdf(0) so should be exactly infinity.
632 negative_binomial_distribution<RealType>(static_cast<RealType>(8), static_cast<RealType>(0.25)),
633 static_cast<RealType>(0))),
634 std::numeric_limits<RealType>::infinity() );
635 } // test for infinity using std::numeric_limits<>::infinity()
636 else
637 { // real_concept case, so check it throws rather than returning infinity.
638 BOOST_CHECK_EQUAL(
639 quantile( // At P == 1 so k failures should be infinite.
640 negative_binomial_distribution<RealType>(static_cast<RealType>(8), static_cast<RealType>(0.25)),
641 static_cast<RealType>(1)),
642 boost::math::tools::max_value<RealType>() );
643
644 BOOST_CHECK_EQUAL(
645 quantile(complement( // Q zero 1 so P == 1 < cdf(0) so should be exactly infinity.
646 negative_binomial_distribution<RealType>(static_cast<RealType>(8), static_cast<RealType>(0.25)),
647 static_cast<RealType>(0))),
648 boost::math::tools::max_value<RealType>());
649 }
650 BOOST_CHECK( // Should work for built-in and real_concept.
651 quantile(complement( // Q very near to 1 so P nearly 1 < so should be large > 384.
652 negative_binomial_distribution<RealType>(static_cast<RealType>(8), static_cast<RealType>(0.25)),
653 static_cast<RealType>(boost::math::tools::min_value<RealType>())))
654 >= static_cast<RealType>(384) );
655
656 BOOST_CHECK_EQUAL(
657 quantile( // P == 0 < cdf(0) so should be zero.
658 negative_binomial_distribution<RealType>(static_cast<RealType>(8), static_cast<RealType>(0.25)),
659 static_cast<RealType>(0)),
660 static_cast<RealType>(0));
661
662 // Quantile Complement boundary cases:
663
664 BOOST_CHECK_EQUAL(
665 quantile(complement( // Q = 1 so P = 0 < cdf(0) so should be exactly zero.
666 negative_binomial_distribution<RealType>(static_cast<RealType>(8), static_cast<RealType>(0.25)),
667 static_cast<RealType>(1))),
668 static_cast<RealType>(0)
669 );
670
671 BOOST_CHECK_EQUAL(
672 quantile(complement( // Q very near 1 so P == epsilon < cdf(0) so should be exactly zero.
673 negative_binomial_distribution<RealType>(static_cast<RealType>(8), static_cast<RealType>(0.25)),
674 static_cast<RealType>(1 - boost::math::tools::epsilon<RealType>()))),
675 static_cast<RealType>(0)
676 );
677
678 // Check that duff arguments throw domain_error:
679 BOOST_CHECK_THROW(
680 pdf( // Negative successes!
681 negative_binomial_distribution<RealType>(static_cast<RealType>(-1), static_cast<RealType>(0.25)),
682 static_cast<RealType>(0)), std::domain_error
683 );
684 BOOST_CHECK_THROW(
685 pdf( // Negative success_fraction!
686 negative_binomial_distribution<RealType>(static_cast<RealType>(8), static_cast<RealType>(-0.25)),
687 static_cast<RealType>(0)), std::domain_error
688 );
689 BOOST_CHECK_THROW(
690 pdf( // Success_fraction > 1!
691 negative_binomial_distribution<RealType>(static_cast<RealType>(8), static_cast<RealType>(1.25)),
692 static_cast<RealType>(0)),
693 std::domain_error
694 );
695 BOOST_CHECK_THROW(
696 pdf( // Negative k argument !
697 negative_binomial_distribution<RealType>(static_cast<RealType>(8), static_cast<RealType>(0.25)),
698 static_cast<RealType>(-1)),
699 std::domain_error
700 );
701 //BOOST_CHECK_THROW(
702 //pdf( // Unlike binomial there is NO limit on k (failures)
703 //negative_binomial_distribution<RealType>(static_cast<RealType>(8), static_cast<RealType>(0.25)),
704 //static_cast<RealType>(9)), std::domain_error
705 //);
706 BOOST_CHECK_THROW(
707 cdf( // Negative k argument !
708 negative_binomial_distribution<RealType>(static_cast<RealType>(8), static_cast<RealType>(0.25)),
709 static_cast<RealType>(-1)),
710 std::domain_error
711 );
712 BOOST_CHECK_THROW(
713 cdf( // Negative success_fraction!
714 negative_binomial_distribution<RealType>(static_cast<RealType>(8), static_cast<RealType>(-0.25)),
715 static_cast<RealType>(0)), std::domain_error
716 );
717 BOOST_CHECK_THROW(
718 cdf( // Success_fraction > 1!
719 negative_binomial_distribution<RealType>(static_cast<RealType>(8), static_cast<RealType>(1.25)),
720 static_cast<RealType>(0)), std::domain_error
721 );
722 BOOST_CHECK_THROW(
723 quantile( // Negative success_fraction!
724 negative_binomial_distribution<RealType>(static_cast<RealType>(8), static_cast<RealType>(-0.25)),
725 static_cast<RealType>(0)), std::domain_error
726 );
727 BOOST_CHECK_THROW(
728 quantile( // Success_fraction > 1!
729 negative_binomial_distribution<RealType>(static_cast<RealType>(8), static_cast<RealType>(1.25)),
730 static_cast<RealType>(0)), std::domain_error
731 );
732 // End of check throwing 'duff' out-of-domain values.
733
734 #define T RealType
735 #include "negative_binomial_quantile.ipp"
736
737 for(unsigned i = 0; i < negative_binomial_quantile_data.size(); ++i)
738 {
739 using namespace boost::math::policies;
740 typedef policy<discrete_quantile<boost::math::policies::real> > P1;
741 typedef policy<discrete_quantile<integer_round_down> > P2;
742 typedef policy<discrete_quantile<integer_round_up> > P3;
743 typedef policy<discrete_quantile<integer_round_outwards> > P4;
744 typedef policy<discrete_quantile<integer_round_inwards> > P5;
745 typedef policy<discrete_quantile<integer_round_nearest> > P6;
746 RealType tol = boost::math::tools::epsilon<RealType>() * 700;
747 if(!boost::is_floating_point<RealType>::value)
748 tol *= 10; // no lanczos approximation implies less accuracy
749 //
750 // Check full real value first:
751 //
752 negative_binomial_distribution<RealType, P1> p1(negative_binomial_quantile_data[i][0], negative_binomial_quantile_data[i][1]);
753 RealType x = quantile(p1, negative_binomial_quantile_data[i][2]);
754 BOOST_CHECK_CLOSE_FRACTION(x, negative_binomial_quantile_data[i][3], tol);
755 x = quantile(complement(p1, negative_binomial_quantile_data[i][2]));
756 BOOST_CHECK_CLOSE_FRACTION(x, negative_binomial_quantile_data[i][4], tol);
757 //
758 // Now with round down to integer:
759 //
760 negative_binomial_distribution<RealType, P2> p2(negative_binomial_quantile_data[i][0], negative_binomial_quantile_data[i][1]);
761 x = quantile(p2, negative_binomial_quantile_data[i][2]);
762 BOOST_CHECK_EQUAL(x, floor(negative_binomial_quantile_data[i][3]));
763 x = quantile(complement(p2, negative_binomial_quantile_data[i][2]));
764 BOOST_CHECK_EQUAL(x, floor(negative_binomial_quantile_data[i][4]));
765 //
766 // Now with round up to integer:
767 //
768 negative_binomial_distribution<RealType, P3> p3(negative_binomial_quantile_data[i][0], negative_binomial_quantile_data[i][1]);
769 x = quantile(p3, negative_binomial_quantile_data[i][2]);
770 BOOST_CHECK_EQUAL(x, ceil(negative_binomial_quantile_data[i][3]));
771 x = quantile(complement(p3, negative_binomial_quantile_data[i][2]));
772 BOOST_CHECK_EQUAL(x, ceil(negative_binomial_quantile_data[i][4]));
773 //
774 // Now with round to integer "outside":
775 //
776 negative_binomial_distribution<RealType, P4> p4(negative_binomial_quantile_data[i][0], negative_binomial_quantile_data[i][1]);
777 x = quantile(p4, negative_binomial_quantile_data[i][2]);
778 BOOST_CHECK_EQUAL(x, negative_binomial_quantile_data[i][2] < 0.5f ? floor(negative_binomial_quantile_data[i][3]) : ceil(negative_binomial_quantile_data[i][3]));
779 x = quantile(complement(p4, negative_binomial_quantile_data[i][2]));
780 BOOST_CHECK_EQUAL(x, negative_binomial_quantile_data[i][2] < 0.5f ? ceil(negative_binomial_quantile_data[i][4]) : floor(negative_binomial_quantile_data[i][4]));
781 //
782 // Now with round to integer "inside":
783 //
784 negative_binomial_distribution<RealType, P5> p5(negative_binomial_quantile_data[i][0], negative_binomial_quantile_data[i][1]);
785 x = quantile(p5, negative_binomial_quantile_data[i][2]);
786 BOOST_CHECK_EQUAL(x, negative_binomial_quantile_data[i][2] < 0.5f ? ceil(negative_binomial_quantile_data[i][3]) : floor(negative_binomial_quantile_data[i][3]));
787 x = quantile(complement(p5, negative_binomial_quantile_data[i][2]));
788 BOOST_CHECK_EQUAL(x, negative_binomial_quantile_data[i][2] < 0.5f ? floor(negative_binomial_quantile_data[i][4]) : ceil(negative_binomial_quantile_data[i][4]));
789 //
790 // Now with round to nearest integer:
791 //
792 negative_binomial_distribution<RealType, P6> p6(negative_binomial_quantile_data[i][0], negative_binomial_quantile_data[i][1]);
793 x = quantile(p6, negative_binomial_quantile_data[i][2]);
794 BOOST_CHECK_EQUAL(x, floor(negative_binomial_quantile_data[i][3] + 0.5f));
795 x = quantile(complement(p6, negative_binomial_quantile_data[i][2]));
796 BOOST_CHECK_EQUAL(x, floor(negative_binomial_quantile_data[i][4] + 0.5f));
797 }
798
799 return;
800 } // template <class RealType> void test_spots(RealType) // Any floating-point type RealType.
801
BOOST_AUTO_TEST_CASE(test_main)802 BOOST_AUTO_TEST_CASE( test_main )
803 {
804 // Check that can generate negative_binomial distribution using the two convenience methods:
805 using namespace boost::math;
806 negative_binomial mynb1(2., 0.5); // Using typedef - default type is double.
807 negative_binomial_distribution<> myf2(2., 0.5); // Using default RealType double.
808
809 // Basic sanity-check spot values.
810
811 // Test some simple double only examples.
812 negative_binomial_distribution<double> my8dist(8., 0.25);
813 // 8 successes (r), 0.25 success fraction = 35% or 1 in 4 successes.
814 // Note: double values (matching the distribution definition) avoid the need for any casting.
815
816 // Check accessor functions return exact values for double at least.
817 BOOST_CHECK_EQUAL(my8dist.successes(), static_cast<double>(8));
818 BOOST_CHECK_EQUAL(my8dist.success_fraction(), static_cast<double>(1./4.));
819
820 // (Parameter value, arbitrarily zero, only communicates the floating point type).
821 #ifdef TEST_FLOAT
822 test_spots(0.0F); // Test float.
823 #endif
824 #ifdef TEST_DOUBLE
825 test_spots(0.0); // Test double.
826 #endif
827 #ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
828 #ifdef TEST_LDOUBLE
829 test_spots(0.0L); // Test long double.
830 #endif
831 #ifndef BOOST_MATH_NO_REAL_CONCEPT_TESTS
832 #ifdef TEST_REAL_CONCEPT
833 test_spots(boost::math::concepts::real_concept(0.)); // Test real concept.
834 #endif
835 #endif
836 #else
837 std::cout << "<note>The long double tests have been disabled on this platform "
838 "either because the long double overloads of the usual math functions are "
839 "not available at all, or because they are too inaccurate for these tests "
840 "to pass.</note>" << std::cout;
841 #endif
842
843
844 } // BOOST_AUTO_TEST_CASE( test_main )
845
846 /*
847
848 Autorun "i:\boost-06-05-03-1300\libs\math\test\Math_test\debug\test_negative_binomial.exe"
849 Running 1 test case...
850 Tolerance = 0.0119209%.
851 Tolerance 5 eps = 5.96046e-007%.
852 Tolerance = 2.22045e-011%.
853 Tolerance 5 eps = 1.11022e-015%.
854 Tolerance = 2.22045e-011%.
855 Tolerance 5 eps = 1.11022e-015%.
856 Tolerance = 2.22045e-011%.
857 Tolerance 5 eps = 1.11022e-015%.
858 *** No errors detected
859
860 */
861