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