]>
Commit | Line | Data |
---|---|---|
7c673cae FG |
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 | |
92f5a8d4 | 41 | #include <boost/test/tools/floating_point_comparison.hpp> // for BOOST_CHECK_CLOSE |
7c673cae FG |
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> | |
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. | |
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 | ||
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 | */ |