]>
git.proxmox.com Git - ceph.git/blob - ceph/src/boost/libs/math/test/test_poisson.cpp
3 // Copyright Paul A. Bristow 2007.
4 // Copyright John Maddock 2006.
6 // Use, modification and distribution are subject to the
7 // Boost Software License, Version 1.0.
8 // (See accompanying file LICENSE_1_0.txt
9 // or copy at http://www.boost.org/LICENSE_1_0.txt)
11 // Basic sanity test for Poisson Cumulative Distribution Function.
13 #define BOOST_MATH_DISCRETE_QUANTILE_POLICY real
15 #if !defined(TEST_FLOAT) && !defined(TEST_DOUBLE) && !defined(TEST_LDOUBLE) && !defined(TEST_REAL_CONCEPT)
19 # define TEST_REAL_CONCEPT
23 # pragma warning(disable: 4127) // conditional expression is constant.
26 #define BOOST_TEST_MAIN
27 #include <boost/test/unit_test.hpp> // Boost.Test
28 #include <boost/test/floating_point_comparison.hpp>
30 #include <boost/math/concepts/real_concept.hpp> // for real_concept
31 #include <boost/math/distributions/poisson.hpp>
32 using boost::math::poisson_distribution
;
33 #include <boost/math/tools/test.hpp> // for real_concept
35 #include <boost/math/special_functions/gamma.hpp> // for (incomplete) gamma.
36 // using boost::math::qamma_Q;
37 #include "table_type.hpp"
38 #include "test_out_of_range.hpp"
43 using std::setprecision
;
47 using std::numeric_limits
;
49 template <class RealType
> // Any floating-point type RealType.
50 void test_spots(RealType
)
52 // Basic sanity checks, tolerance is about numeric_limits<RealType>::digits10 decimal places,
53 // guaranteed for type RealType, eg 6 for float, 15 for double,
54 // expressed as a percentage (so -2) for BOOST_CHECK_CLOSE,
56 int decdigits
= numeric_limits
<RealType
>::digits10
;
57 // May eb >15 for 80 and 128-bit FP typtes.
59 { // decdigits is not defined, for example real concept,
60 // so assume precision of most test data is double (for example, MathCAD).
61 decdigits
= numeric_limits
<double>::digits10
; // == 15 for 64-bit
63 if (decdigits
> 15 ) // numeric_limits<double>::digits10)
64 { // 15 is the accuracy of the MathCAD test data.
65 decdigits
= 15; // numeric_limits<double>::digits10;
68 decdigits
-= 1; // Perhaps allow some decimal digit(s) margin of numerical error.
69 RealType tolerance
= static_cast<RealType
>(std::pow(10., static_cast<double>(2-decdigits
))); // 1e-6 (-2 so as %)
70 tolerance
*= 2; // Allow some bit(s) small margin (2 means + or - 1 bit) of numerical error.
71 // Typically 2e-13% = 2e-15 as fraction for double.
73 // Sources of spot test values:
75 // Many be some combinations for which the result is 'exact',
76 // or at least is good to 40 decimal digits.
77 // 40 decimal digits includes 128-bit significand User Defined Floating-Point types,
79 // Best source of accurate values is:
80 // Mathworld online calculator (40 decimal digits precision, suitable for up to 128-bit significands)
81 // http://functions.wolfram.com/webMathematica/FunctionEvaluation.jsp?name=GammaRegularized
82 // GammaRegularized is same as gamma incomplete, gamma or gamma_q(a, x) or Q(a, z).
84 // http://documents.wolfram.com/calculationcenter/v2/Functions/ListsMatrices/Statistics/PoissonDistribution.html
86 // MathCAD defines ppois(k, lambda== mean) as k integer, k >=0.
87 // ppois(0, 5) = 6.73794699908547e-3
88 // ppois(1, 5) = 0.040427681994513;
89 // ppois(10, 10) = 5.830397501929850E-001
90 // ppois(10, 1) = 9.999999899522340E-001
91 // ppois(5,5) = 0.615960654833065
93 // qpois returns inverse Poission distribution, that is the smallest (floor) k so that ppois(k, lambda) >= p
94 // p is real number, real mean lambda > 0
95 // k is approximately the integer for which probability(X <= k) = p
96 // when random variable X has the Poisson distribution with parameters lambda.
97 // Uses discrete bisection.
98 // qpois(6.73794699908547e-3, 5) = 1
99 // qpois(0.040427681994513, 5) =
101 // Test Poisson with spot values from MathCAD 'known good'.
103 using boost::math::poisson_distribution
;
104 using ::boost::math::poisson
;
105 using ::boost::math::cdf
;
106 using ::boost::math::pdf
;
108 // Check that bad arguments throw.
109 BOOST_MATH_CHECK_THROW(
110 cdf(poisson_distribution
<RealType
>(static_cast<RealType
>(0)), // mean zero is bad.
111 static_cast<RealType
>(0)), // even for a good k.
112 std::domain_error
); // Expected error to be thrown.
114 BOOST_MATH_CHECK_THROW(
115 cdf(poisson_distribution
<RealType
>(static_cast<RealType
>(-1)), // mean negative is bad.
116 static_cast<RealType
>(0)),
119 BOOST_MATH_CHECK_THROW(
120 cdf(poisson_distribution
<RealType
>(static_cast<RealType
>(1)), // mean unit OK,
121 static_cast<RealType
>(-1)), // but negative events is bad.
124 BOOST_MATH_CHECK_THROW(
125 cdf(poisson_distribution
<RealType
>(static_cast<RealType
>(0)), // mean zero is bad.
126 static_cast<RealType
>(99999)), // for any k events.
129 BOOST_MATH_CHECK_THROW(
130 cdf(poisson_distribution
<RealType
>(static_cast<RealType
>(0)), // mean zero is bad.
131 static_cast<RealType
>(99999)), // for any k events.
134 BOOST_MATH_CHECK_THROW(
135 quantile(poisson_distribution
<RealType
>(static_cast<RealType
>(0)), // mean zero.
136 static_cast<RealType
>(0.5)), // probability OK.
139 BOOST_MATH_CHECK_THROW(
140 quantile(poisson_distribution
<RealType
>(static_cast<RealType
>(-1)),
141 static_cast<RealType
>(-1)), // bad probability.
144 BOOST_MATH_CHECK_THROW(
145 quantile(poisson_distribution
<RealType
>(static_cast<RealType
>(1)),
146 static_cast<RealType
>(-1)), // bad probability.
149 BOOST_MATH_CHECK_THROW(
150 quantile(poisson_distribution
<RealType
>(static_cast<RealType
>(1)),
151 static_cast<RealType
>(1)), // bad probability.
152 std::overflow_error
);
154 BOOST_MATH_CHECK_THROW(
155 quantile(complement(poisson_distribution
<RealType
>(static_cast<RealType
>(1)),
156 static_cast<RealType
>(0))), // bad probability.
157 std::overflow_error
);
160 quantile(poisson_distribution
<RealType
>(static_cast<RealType
>(1)),
161 static_cast<RealType
>(0)), // bad probability.
165 quantile(complement(poisson_distribution
<RealType
>(static_cast<RealType
>(1)),
166 static_cast<RealType
>(1))), // bad probability.
169 // Check some test values.
171 BOOST_CHECK_CLOSE( // mode
172 mode(poisson_distribution
<RealType
>(static_cast<RealType
>(4))), // mode = mean = 4.
173 static_cast<RealType
>(4), // mode.
176 //BOOST_CHECK_CLOSE( // mode
177 // median(poisson_distribution<RealType>(static_cast<RealType>(4))), // mode = mean = 4.
178 // static_cast<RealType>(4), // mode.
180 poisson_distribution
<RealType
> dist4(static_cast<RealType
>(40));
182 BOOST_CHECK_CLOSE( // median
183 median(dist4
), // mode = mean = 4. median = 40.328333333333333
184 quantile(dist4
, static_cast<RealType
>(0.5)), // 39.332839138842637
189 pdf(poisson_distribution
<RealType
>(static_cast<RealType
>(4)), // mean 4.
190 static_cast<RealType
>(0)),
191 static_cast<RealType
>(1.831563888873410E-002), // probability.
195 pdf(poisson_distribution
<RealType
>(static_cast<RealType
>(4)), // mean 4.
196 static_cast<RealType
>(2)),
197 static_cast<RealType
>(1.465251111098740E-001), // probability.
201 pdf(poisson_distribution
<RealType
>(static_cast<RealType
>(20)), // mean big.
202 static_cast<RealType
>(1)), // k small
203 static_cast<RealType
>(4.122307244877130E-008), // probability.
207 pdf(poisson_distribution
<RealType
>(static_cast<RealType
>(4)), // mean 4.
208 static_cast<RealType
>(20)), // K>> mean
209 static_cast<RealType
>(8.277463646553730E-009), // probability.
214 cdf(poisson_distribution
<RealType
>(static_cast<RealType
>(1)), // mean unity.
215 static_cast<RealType
>(0)), // zero k events.
216 static_cast<RealType
>(3.678794411714420E-1), // probability.
220 cdf(poisson_distribution
<RealType
>(static_cast<RealType
>(1)), // mean unity.
221 static_cast<RealType
>(1)), // one k event.
222 static_cast<RealType
>(7.357588823428830E-1), // probability.
226 cdf(poisson_distribution
<RealType
>(static_cast<RealType
>(1)), // mean unity.
227 static_cast<RealType
>(2)), // two k events.
228 static_cast<RealType
>(9.196986029286060E-1), // probability.
232 cdf(poisson_distribution
<RealType
>(static_cast<RealType
>(1)), // mean unity.
233 static_cast<RealType
>(10)), // two k events.
234 static_cast<RealType
>(9.999999899522340E-1), // probability.
238 cdf(poisson_distribution
<RealType
>(static_cast<RealType
>(1)), // mean unity.
239 static_cast<RealType
>(15)), // two k events.
240 static_cast<RealType
>(9.999999999999810E-1), // probability.
244 cdf(poisson_distribution
<RealType
>(static_cast<RealType
>(1)), // mean unity.
245 static_cast<RealType
>(16)), // two k events.
246 static_cast<RealType
>(9.999999999999990E-1), // probability.
250 cdf(poisson_distribution
<RealType
>(static_cast<RealType
>(1)), // mean unity.
251 static_cast<RealType
>(17)), // two k events.
252 static_cast<RealType
>(1.), // probability unity for double.
256 cdf(poisson_distribution
<RealType
>(static_cast<RealType
>(1)), // mean unity.
257 static_cast<RealType
>(33)), // k events at limit for float unchecked_factorial table.
258 static_cast<RealType
>(1.), // probability.
262 cdf(poisson_distribution
<RealType
>(static_cast<RealType
>(100)), // mean 100.
263 static_cast<RealType
>(33)), // k events at limit for float unchecked_factorial table.
264 static_cast<RealType
>(6.328271240363390E-15), // probability is tiny.
265 tolerance
* static_cast<RealType
>(2e11
)); // 6.3495253382825722e-015 MathCAD
266 // Note that there two tiny probability are much more different.
269 cdf(poisson_distribution
<RealType
>(static_cast<RealType
>(100)), // mean 100.
270 static_cast<RealType
>(34)), // k events at limit for float unchecked_factorial table.
271 static_cast<RealType
>(1.898481372109020E-14), // probability is tiny.
272 tolerance
*static_cast<RealType
>(2e11
)); // 1.8984813721090199e-014 MathCAD
276 cdf(poisson_distribution
<RealType
>(static_cast<RealType
>(33)), // mean = k
277 static_cast<RealType
>(33)), // k events above limit for float unchecked_factorial table.
278 static_cast<RealType
>(5.461191812386560E-1), // probability.
282 cdf(poisson_distribution
<RealType
>(static_cast<RealType
>(33)), // mean = k-1
283 static_cast<RealType
>(34)), // k events above limit for float unchecked_factorial table.
284 static_cast<RealType
>(6.133535681502950E-1), // probability.
288 cdf(poisson_distribution
<RealType
>(static_cast<RealType
>(1)), // mean unity.
289 static_cast<RealType
>(34)), // k events above limit for float unchecked_factorial table.
290 static_cast<RealType
>(1.), // probability.
294 cdf(poisson_distribution
<RealType
>(static_cast<RealType
>(5.)), // mean
295 static_cast<RealType
>(5)), // k events.
296 static_cast<RealType
>(0.615960654833065), // probability.
299 cdf(poisson_distribution
<RealType
>(static_cast<RealType
>(5.)), // mean
300 static_cast<RealType
>(1)), // k events.
301 static_cast<RealType
>(0.040427681994512805), // probability.
305 cdf(poisson_distribution
<RealType
>(static_cast<RealType
>(5.)), // mean
306 static_cast<RealType
>(0)), // k events (uses special case formula, not gamma).
307 static_cast<RealType
>(0.006737946999085467), // probability.
311 cdf(poisson_distribution
<RealType
>(static_cast<RealType
>(1.)), // mean
312 static_cast<RealType
>(0)), // k events (uses special case formula, not gamma).
313 static_cast<RealType
>(0.36787944117144233), // probability.
317 cdf(poisson_distribution
<RealType
>(static_cast<RealType
>(10.)), // mean
318 static_cast<RealType
>(10)), // k events.
319 static_cast<RealType
>(0.5830397501929856), // probability.
323 cdf(poisson_distribution
<RealType
>(static_cast<RealType
>(4.)), // mean
324 static_cast<RealType
>(5)), // k events.
325 static_cast<RealType
>(0.785130387030406), // probability.
329 BOOST_CHECK_CLOSE( // Complement CDF
330 cdf(complement(poisson_distribution
<RealType
>(static_cast<RealType
>(4.)), // mean
331 static_cast<RealType
>(5))), // k events.
332 static_cast<RealType
>(1 - 0.785130387030406), // probability.
335 BOOST_CHECK_CLOSE( // Complement CDF
336 cdf(complement(poisson_distribution
<RealType
>(static_cast<RealType
>(4.)), // mean
337 static_cast<RealType
>(0))), // Zero k events (uses special case formula, not gamma).
338 static_cast<RealType
>(0.98168436111126578), // probability.
340 BOOST_CHECK_CLOSE( // Complement CDF
341 cdf(complement(poisson_distribution
<RealType
>(static_cast<RealType
>(1.)), // mean
342 static_cast<RealType
>(0))), // Zero k events (uses special case formula, not gamma).
343 static_cast<RealType
>(0.63212055882855767), // probability.
346 // Example where k is bigger than max_factorial (>34 for float)
347 // (therefore using log gamma so perhaps less accurate).
349 cdf(poisson_distribution
<RealType
>(static_cast<RealType
>(40.)), // mean
350 static_cast<RealType
>(40)), // k events.
351 static_cast<RealType
>(0.5419181783625430), // probability.
354 // Quantile & complement.
356 boost::math::quantile(
357 poisson_distribution
<RealType
>(5), // mean.
358 static_cast<RealType
>(0.615960654833065)), // probability.
359 static_cast<RealType
>(5.), // Expect k = 5
362 // EQUAL is too optimistic - fails [5.0000000000000124 != 5]
363 // BOOST_CHECK_EQUAL(boost::math::quantile( //
364 // poisson_distribution<RealType>(5.), // mean.
365 // static_cast<RealType>(0.615960654833065)), // probability.
366 // static_cast<RealType>(5.)); // Expect k = 5 events.
368 BOOST_CHECK_CLOSE(boost::math::quantile(
369 poisson_distribution
<RealType
>(4), // mean.
370 static_cast<RealType
>(0.785130387030406)), // probability.
371 static_cast<RealType
>(5.), // Expect k = 5 events.
374 // Check on quantile of other examples of inverse of cdf.
376 cdf(poisson_distribution
<RealType
>(static_cast<RealType
>(10.)), // mean
377 static_cast<RealType
>(10)), // k events.
378 static_cast<RealType
>(0.5830397501929856), // probability.
381 BOOST_CHECK_CLOSE(boost::math::quantile( // inverse of cdf above.
382 poisson_distribution
<RealType
>(10.), // mean.
383 static_cast<RealType
>(0.5830397501929856)), // probability.
384 static_cast<RealType
>(10.), // Expect k = 10 events.
389 cdf(poisson_distribution
<RealType
>(static_cast<RealType
>(4.)), // mean
390 static_cast<RealType
>(5)), // k events.
391 static_cast<RealType
>(0.785130387030406), // probability.
394 BOOST_CHECK_CLOSE(boost::math::quantile( // inverse of cdf above.
395 poisson_distribution
<RealType
>(4.), // mean.
396 static_cast<RealType
>(0.785130387030406)), // probability.
397 static_cast<RealType
>(5.), // Expect k = 10 events.
402 //BOOST_CHECK_CLOSE(boost::math::quantile(
403 // poisson_distribution<RealType>(5), // mean.
404 // static_cast<RealType>(0.785130387030406)), // probability.
405 // // 6.1882832344329559 result but MathCAD givest smallest integer ppois(k, mean) >= prob
406 // static_cast<RealType>(6.), // Expect k = 6 events.
409 //BOOST_CHECK_CLOSE(boost::math::quantile(
410 // poisson_distribution<RealType>(5), // mean.
411 // static_cast<RealType>(0.77)), // probability.
412 // // 6.1882832344329559 result but MathCAD givest smallest integer ppois(k, mean) >= prob
413 // static_cast<RealType>(7.), // Expect k = 6 events.
416 //BOOST_CHECK_CLOSE(boost::math::quantile(
417 // poisson_distribution<RealType>(5), // mean.
418 // static_cast<RealType>(0.75)), // probability.
419 // // 6.1882832344329559 result but MathCAD givest smallest integer ppois(k, mean) >= prob
420 // static_cast<RealType>(6.), // Expect k = 6 events.
424 boost::math::quantile(
426 poisson_distribution
<RealType
>(4),
427 static_cast<RealType
>(1 - 0.785130387030406))), // complement.
428 static_cast<RealType
>(5), // Expect k = 5 events.
431 BOOST_CHECK_EQUAL(boost::math::quantile( // Check case when probability < cdf(0) (== pdf(0))
432 poisson_distribution
<RealType
>(1), // mean is small, so cdf and pdf(0) are about 0.35.
433 static_cast<RealType
>(0.0001)), // probability < cdf(0).
434 static_cast<RealType
>(0)); // Expect k = 0 events exactly.
437 boost::math::quantile(
439 poisson_distribution
<RealType
>(1),
440 static_cast<RealType
>(0.9999))), // complement, so 1-probability < cdf(0)
441 static_cast<RealType
>(0)); // Expect k = 0 events exactly.
444 // Test quantile policies against test data:
447 #include "poisson_quantile.ipp"
449 for(unsigned i
= 0; i
< poisson_quantile_data
.size(); ++i
)
451 using namespace boost::math::policies
;
452 typedef policy
<discrete_quantile
<real
> > P1
;
453 typedef policy
<discrete_quantile
<integer_round_down
> > P2
;
454 typedef policy
<discrete_quantile
<integer_round_up
> > P3
;
455 typedef policy
<discrete_quantile
<integer_round_outwards
> > P4
;
456 typedef policy
<discrete_quantile
<integer_round_inwards
> > P5
;
457 typedef policy
<discrete_quantile
<integer_round_nearest
> > P6
;
458 RealType tol
= boost::math::tools::epsilon
<RealType
>() * 20;
459 if(!boost::is_floating_point
<RealType
>::value
)
462 // Check full real value first:
464 poisson_distribution
<RealType
, P1
> p1(poisson_quantile_data
[i
][0]);
465 RealType x
= quantile(p1
, poisson_quantile_data
[i
][1]);
466 BOOST_CHECK_CLOSE_FRACTION(x
, poisson_quantile_data
[i
][2], tol
);
467 x
= quantile(complement(p1
, poisson_quantile_data
[i
][1]));
468 BOOST_CHECK_CLOSE_FRACTION(x
, poisson_quantile_data
[i
][3], tol
* 3);
470 // Now with round down to integer:
472 poisson_distribution
<RealType
, P2
> p2(poisson_quantile_data
[i
][0]);
473 x
= quantile(p2
, poisson_quantile_data
[i
][1]);
474 BOOST_CHECK_EQUAL(x
, floor(poisson_quantile_data
[i
][2]));
475 x
= quantile(complement(p2
, poisson_quantile_data
[i
][1]));
476 BOOST_CHECK_EQUAL(x
, floor(poisson_quantile_data
[i
][3]));
478 // Now with round up to integer:
480 poisson_distribution
<RealType
, P3
> p3(poisson_quantile_data
[i
][0]);
481 x
= quantile(p3
, poisson_quantile_data
[i
][1]);
482 BOOST_CHECK_EQUAL(x
, ceil(poisson_quantile_data
[i
][2]));
483 x
= quantile(complement(p3
, poisson_quantile_data
[i
][1]));
484 BOOST_CHECK_EQUAL(x
, ceil(poisson_quantile_data
[i
][3]));
486 // Now with round to integer "outside":
488 poisson_distribution
<RealType
, P4
> p4(poisson_quantile_data
[i
][0]);
489 x
= quantile(p4
, poisson_quantile_data
[i
][1]);
490 BOOST_CHECK_EQUAL(x
, poisson_quantile_data
[i
][1] < 0.5f
? floor(poisson_quantile_data
[i
][2]) : ceil(poisson_quantile_data
[i
][2]));
491 x
= quantile(complement(p4
, poisson_quantile_data
[i
][1]));
492 BOOST_CHECK_EQUAL(x
, poisson_quantile_data
[i
][1] < 0.5f
? ceil(poisson_quantile_data
[i
][3]) : floor(poisson_quantile_data
[i
][3]));
494 // Now with round to integer "inside":
496 poisson_distribution
<RealType
, P5
> p5(poisson_quantile_data
[i
][0]);
497 x
= quantile(p5
, poisson_quantile_data
[i
][1]);
498 BOOST_CHECK_EQUAL(x
, poisson_quantile_data
[i
][1] < 0.5f
? ceil(poisson_quantile_data
[i
][2]) : floor(poisson_quantile_data
[i
][2]));
499 x
= quantile(complement(p5
, poisson_quantile_data
[i
][1]));
500 BOOST_CHECK_EQUAL(x
, poisson_quantile_data
[i
][1] < 0.5f
? floor(poisson_quantile_data
[i
][3]) : ceil(poisson_quantile_data
[i
][3]));
502 // Now with round to nearest integer:
504 poisson_distribution
<RealType
, P6
> p6(poisson_quantile_data
[i
][0]);
505 x
= quantile(p6
, poisson_quantile_data
[i
][1]);
506 BOOST_CHECK_EQUAL(x
, floor(poisson_quantile_data
[i
][2] + 0.5f
));
507 x
= quantile(complement(p6
, poisson_quantile_data
[i
][1]));
508 BOOST_CHECK_EQUAL(x
, floor(poisson_quantile_data
[i
][3] + 0.5f
));
510 check_out_of_range
<poisson_distribution
<RealType
> >(1);
511 } // template <class RealType>void test_spots(RealType)
515 BOOST_AUTO_TEST_CASE( test_main
)
517 // Check that can construct normal distribution using the two convenience methods:
518 using namespace boost::math
;
519 poisson
myp1(2); // Using typedef
520 poisson_distribution
<> myp2(2); // Using default RealType double.
522 // Basic sanity-check spot values.
524 // Some plain double examples & tests:
525 cout
.precision(17); // double max_digits10
526 cout
.setf(ios::showpoint
);
528 poisson
mypoisson(4.); // // mean = 4, default FP type is double.
529 cout
<< "mean(mypoisson, 4.) == " << mean(mypoisson
) << endl
;
530 cout
<< "mean(mypoisson, 0.) == " << mean(mypoisson
) << endl
;
531 cout
<< "cdf(mypoisson, 2.) == " << cdf(mypoisson
, 2.) << endl
;
532 cout
<< "pdf(mypoisson, 2.) == " << pdf(mypoisson
, 2.) << endl
;
534 // poisson mydudpoisson(0.);
535 // throws (if BOOST_MATH_DOMAIN_ERROR_POLICY == throw_on_error).
537 #ifndef BOOST_NO_EXCEPTIONS
538 BOOST_MATH_CHECK_THROW(poisson
mydudpoisson(-1), std::domain_error
);// Mean must be > 0.
539 BOOST_MATH_CHECK_THROW(poisson
mydudpoisson(-1), std::logic_error
);// Mean must be > 0.
541 BOOST_MATH_CHECK_THROW(poisson(-1), std::domain_error
);// Mean must be > 0.
542 BOOST_MATH_CHECK_THROW(poisson(-1), std::logic_error
);// Mean must be > 0.
544 // Passes the check because logic_error is a parent????
545 // BOOST_MATH_CHECK_THROW(poisson mydudpoisson(-1), std::overflow_error); // fails the check
546 // because overflow_error is unrelated - except from std::exception
547 BOOST_MATH_CHECK_THROW(cdf(mypoisson
, -1), std::domain_error
); // k must be >= 0
549 BOOST_CHECK_EQUAL(mean(mypoisson
), 4.);
551 pdf(mypoisson
, 2.), // k events = 2.
552 1.465251111098740E-001, // probability.
556 cdf(mypoisson
, 2.), // k events = 2.
557 0.238103305553545, // probability.
562 // Compare cdf from finite sum of pdf and gamma_q.
563 using boost::math::cdf
;
564 using boost::math::pdf
;
567 cout
.precision(17); // double max_digits10
568 cout
.setf(ios::showpoint
);
569 cout
<< showpoint
<< endl
; // Ensure trailing zeros are shown.
570 // This also helps show the expected precision max_digits10
571 //cout.unsetf(ios::showpoint); // No trailing zeros are shown.
573 cout
<< "k pdf sum cdf diff" << endl
;
575 for (int i
= 0; i
<= 50; i
++)
578 double p
= pdf(poisson_distribution
<double>(mean
), static_cast<double>(i
));
581 cout
<< p
<< ' ' << sum
<< ' '
582 << cdf(poisson_distribution
<double>(mean
), static_cast<double>(i
)) << ' ';
584 cout
<< boost::math::gamma_q
<double>(i
+1, mean
); // cdf
585 double diff
= boost::math::gamma_q
<double>(i
+1, mean
) - sum
; // cdf -sum
586 cout
<< setprecision (2) << ' ' << diff
; // 0 0 to 4, 1 eps 5 to 9, 10 to 20 2 eps, 21 upwards 3 eps
590 cdf(mypoisson
, static_cast<double>(i
)),
592 4e-14); // Fails at 2e-14
593 // This call puts the precision etc back to default 6 !!!
594 cout
<< setprecision(17) << showpoint
;
600 cout
<< cdf(poisson_distribution
<double>(5), static_cast<double>(0)) << ' ' << endl
; // 0.006737946999085467
601 cout
<< cdf(poisson_distribution
<double>(5), static_cast<double>(1)) << ' ' << endl
; // 0.040427681994512805
602 cout
<< cdf(poisson_distribution
<double>(2), static_cast<double>(3)) << ' ' << endl
; // 0.85712346049854715
604 { // Compare approximate formula in Wikipedia with quantile(half)
605 for (int i
= 1; i
< 100; i
++)
607 poisson_distribution
<double> distn(static_cast<double>(i
));
608 cout
<< i
<< ' ' << median(distn
) << ' ' << quantile(distn
, 0.5) << ' '
609 << median(distn
) - quantile(distn
, 0.5) << endl
; // formula appears to be out-by-one??
610 } // so quantile(half) used via derived accressors.
614 // (Parameter value, arbitrarily zero, only communicates the floating-point type).
616 test_spots(0.0F
); // Test float.
619 test_spots(0.0); // Test double.
621 #ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
622 if (numeric_limits
<long double>::digits10
> numeric_limits
<double>::digits10
)
623 { // long double is better than double (so not MSVC where they are same).
625 test_spots(0.0L); // Test long double.
629 #ifndef BOOST_MATH_NO_REAL_CONCEPT_TESTS
630 #ifdef TEST_REAL_CONCEPT
631 test_spots(boost::math::concepts::real_concept(0.)); // Test real concept.
636 } // BOOST_AUTO_TEST_CASE( test_main )
642 Autorun "i:\boost-06-05-03-1300\libs\math\test\Math_test\debug\test_poisson.exe"
643 Running 1 test case...
644 mean(mypoisson, 4.) == 4.0000000000000000
645 mean(mypoisson, 0.) == 4.0000000000000000
646 cdf(mypoisson, 2.) == 0.23810330555354431
647 pdf(mypoisson, 2.) == 0.14652511110987343
648 *** No errors detected