1 // test_geometric.cpp
2
3 // Copyright Paul A. Bristow 2010.
4 // Copyright John Maddock 2010.
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 Geometric Distribution.
12
13 // Note that these defines must be placed BEFORE #includes.
14 #define BOOST_MATH_OVERFLOW_ERROR_POLICY ignore_error
15 // because several tests overflow & underflow by design.
16 #define BOOST_MATH_DISCRETE_QUANTILE_POLICY real
17
18 #ifdef _MSC_VER
19 # pragma warning(disable: 4127) // conditional expression is constant.
20 #endif
21
22 #if !defined(TEST_FLOAT) && !defined(TEST_DOUBLE) && !defined(TEST_LDOUBLE) && !defined(TEST_REAL_CONCEPT)
23 # define TEST_FLOAT
24 # define TEST_DOUBLE
25 # define TEST_LDOUBLE
26 # define TEST_REAL_CONCEPT
27 #endif
28
29 #include <boost/math/concepts/real_concept.hpp> // for real_concept
30 using ::boost::math::concepts::real_concept;
31
32 #include <boost/math/distributions/geometric.hpp> // for geometric_distribution
33 using boost::math::geometric_distribution;
34 using boost::math::geometric; // using typedef for geometric_distribution<double>
35
36 #include <boost/math/distributions/negative_binomial.hpp> // for some comparisons.
37
38 #define BOOST_TEST_MAIN
39 #include <boost/test/unit_test.hpp> // for test_main
40 #include <boost/test/floating_point_comparison.hpp> // for BOOST_CHECK_CLOSE_FRACTION
41 #include "test_out_of_range.hpp"
42
43 #include <iostream>
44 using std::cout;
45 using std::endl;
46 using std::setprecision;
47 using std::showpoint;
48 #include <limits>
49 using std::numeric_limits;
50
51 template <class RealType>
test_spot(RealType k,RealType p,RealType P,RealType Q,RealType tol)52 void test_spot( // Test a single spot value against 'known good' values.
53 RealType k, // Number of failures.
54 RealType p, // Probability of success_fraction.
55 RealType P, // CDF probability.
56 RealType Q, // Complement of CDF.
57 RealType tol) // Test tolerance.
58 {
59 boost::math::geometric_distribution<RealType> g(p);
60 BOOST_CHECK_EQUAL(p, g.success_fraction());
61 BOOST_CHECK_CLOSE_FRACTION(cdf(g, k), P, tol);
62
63 if((P < 0.99) && (Q < 0.99))
64 {
65 // We can only check this if P is not too close to 1,
66 // so that we can guarantee that Q is free of error:
67 //
68 BOOST_CHECK_CLOSE_FRACTION(
69 cdf(complement(g, k)), Q, tol);
70 if(k != 0)
71 {
72 BOOST_CHECK_CLOSE_FRACTION(
73 quantile(g, P), k, tol);
74 }
75 else
76 {
77 // Just check quantile is very small:
78 if((std::numeric_limits<RealType>::max_exponent <= std::numeric_limits<double>::max_exponent)
79 && (boost::is_floating_point<RealType>::value))
80 {
81 // Limit where this is checked: if exponent range is very large we may
82 // run out of iterations in our root finding algorithm.
83 BOOST_CHECK(quantile(g, P) < boost::math::tools::epsilon<RealType>() * 10);
84 }
85 }
86 if(k != 0)
87 {
88 BOOST_CHECK_CLOSE_FRACTION(
89 quantile(complement(g, Q)), k, tol);
90 }
91 else
92 {
93 // Just check quantile is very small:
94 if((std::numeric_limits<RealType>::max_exponent <= std::numeric_limits<double>::max_exponent)
95 && (boost::is_floating_point<RealType>::value))
96 {
97 // Limit where this is checked: if exponent range is very large we may
98 // run out of iterations in our root finding algorithm.
99 BOOST_CHECK(quantile(complement(g, Q)) < boost::math::tools::epsilon<RealType>() * 10);
100 }
101 }
102 } // if((P < 0.99) && (Q < 0.99))
103
104 // Parameter estimation test: estimate success ratio:
105 BOOST_CHECK_CLOSE_FRACTION(
106 geometric_distribution<RealType>::find_lower_bound_on_p(
107 1+k, P),
108 p, 0.02); // Wide tolerance needed for some tests.
109 // Note we bump up the sample size here, purely for the sake of the test,
110 // internally the function has to adjust the sample size so that we get
111 // the right upper bound, our test undoes this, so we can verify the result.
112 BOOST_CHECK_CLOSE_FRACTION(
113 geometric_distribution<RealType>::find_upper_bound_on_p(
114 1+k+1, Q),
115 p, 0.02);
116
117 if(Q < P)
118 {
119 //
120 // We check two things here, that the upper and lower bounds
121 // are the right way around, and that they do actually bracket
122 // the naive estimate of p = successes / (sample size)
123 //
124 BOOST_CHECK(
125 geometric_distribution<RealType>::find_lower_bound_on_p(
126 1+k, Q)
127 <=
128 geometric_distribution<RealType>::find_upper_bound_on_p(
129 1+k, Q)
130 );
131 BOOST_CHECK(
132 geometric_distribution<RealType>::find_lower_bound_on_p(
133 1+k, Q)
134 <=
135 1 / (1+k)
136 );
137 BOOST_CHECK(
138 1 / (1+k)
139 <=
140 geometric_distribution<RealType>::find_upper_bound_on_p(
141 1+k, Q)
142 );
143 }
144 else
145 {
146 // As above but when P is small.
147 BOOST_CHECK(
148 geometric_distribution<RealType>::find_lower_bound_on_p(
149 1+k, P)
150 <=
151 geometric_distribution<RealType>::find_upper_bound_on_p(
152 1+k, P)
153 );
154 BOOST_CHECK(
155 geometric_distribution<RealType>::find_lower_bound_on_p(
156 1+k, P)
157 <=
158 1 / (1+k)
159 );
160 BOOST_CHECK(
161 1 / (1+k)
162 <=
163 geometric_distribution<RealType>::find_upper_bound_on_p(
164 1+k, P)
165 );
166 }
167
168 // Estimate sample size:
169 BOOST_CHECK_CLOSE_FRACTION(
170 geometric_distribution<RealType>::find_minimum_number_of_trials(
171 k, p, P),
172 1+k, 0.02); // Can differ 50 to 51 for small p
173 BOOST_CHECK_CLOSE_FRACTION(
174 geometric_distribution<RealType>::find_maximum_number_of_trials(
175 k, p, Q),
176 1+k, 0.02);
177
178 } // test_spot
179
180 template <class RealType> // Any floating-point type RealType.
test_spots(RealType)181 void test_spots(RealType)
182 {
183 // Basic sanity checks.
184 // Most test data is to double precision (17 decimal digits) only,
185
186 cout << "Floating point Type is " << typeid(RealType).name() << endl;
187
188 // so set tolerance to 1000 eps expressed as a fraction,
189 // or 1000 eps of type double expressed as a fraction,
190 // whichever is the larger.
191
192 RealType tolerance = (std::max)
193 (boost::math::tools::epsilon<RealType>(),
194 static_cast<RealType>(std::numeric_limits<double>::epsilon()));
195 tolerance *= 10; // 10 eps
196
197 cout << "Tolerance = " << tolerance << "." << endl;
198
199 RealType tol1eps = boost::math::tools::epsilon<RealType>(); // Very tight, suit exact values.
200 //RealType tol2eps = boost::math::tools::epsilon<RealType>() * 2; // Tight, values.
201 RealType tol5eps = boost::math::tools::epsilon<RealType>() * 5; // Wider 5 epsilon.
202 cout << "Tolerance 5 eps = " << tol5eps << "." << endl;
203
204
205 // Sources of spot test values are mainly R.
206
207 using boost::math::geometric_distribution;
208 using boost::math::geometric;
209 using boost::math::cdf;
210 using boost::math::pdf;
211 using boost::math::quantile;
212 using boost::math::complement;
213
214 BOOST_MATH_STD_USING // for std math functions
215
216 // Test geometric using cdf spot values R
217 // These test quantiles and complements as well.
218
219 test_spot( //
220 static_cast<RealType>(2), // Number of failures, k
221 static_cast<RealType>(0.5), // Probability of success as fraction, p
222 static_cast<RealType>(0.875L), // Probability of result (CDF), P
223 static_cast<RealType>(0.125L), // complement CCDF Q = 1 - P
224 tolerance);
225
226 test_spot( //
227 static_cast<RealType>(0), // Number of failures, k
228 static_cast<RealType>(0.25), // Probability of success as fraction, p
229 static_cast<RealType>(0.25), // Probability of result (CDF), P
230 static_cast<RealType>(0.75), // Q = 1 - P
231 tolerance);
232
233 test_spot(
234 // R formatC(pgeom(10,0.25), digits=17) [1] "0.95776486396789551"
235 // formatC(pgeom(10,0.25, FALSE), digits=17) [1] "0.042235136032104499"
236
237 static_cast<RealType>(10), // Number of failures, k
238 static_cast<RealType>(0.25), // Probability of success, p
239 static_cast<RealType>(0.95776486396789551L), // Probability of result (CDF), P
240 static_cast<RealType>(0.042235136032104499L), // Q = 1 - P
241 tolerance);
242
243 test_spot( //
244 // > R formatC(pgeom(50,0.25, TRUE), digits=17) [1] "0.99999957525875771"
245 // > R formatC(pgeom(50,0.25, FALSE), digits=17) [1] "4.2474124232020353e-07"
246 static_cast<RealType>(50), // Number of failures, k
247 static_cast<RealType>(0.25), // Probability of success, p
248 static_cast<RealType>(0.99999957525875771), // Probability of result (CDF), P
249 static_cast<RealType>(4.2474124232020353e-07), // Q = 1 - P
250 tolerance);
251 /*
252 // This causes failures in find_upper_bound_on_p p is small branch.
253 test_spot( // formatC(pgeom(50,0.01, TRUE), digits=17)[1] "0.40104399353383874"
254 // > formatC(pgeom(50,0.01, FALSE), digits=17) [1] "0.59895600646616121"
255 static_cast<RealType>(50), // Number of failures, k
256 static_cast<RealType>(0.01), // Probability of success, p
257 static_cast<RealType>(0.40104399353383874), // Probability of result (CDF), P
258 static_cast<RealType>(0.59895600646616121), // Q = 1 - P
259 tolerance);
260 */
261
262 test_spot( // > formatC(pgeom(50,0.99, TRUE), digits=17) [1] " 1"
263 // formatC(pgeom(50,0.99, FALSE), digits=17) [1] "1.0000000000000364e-102"
264 static_cast<RealType>(50), // Number of failures, k
265 static_cast<RealType>(0.99), // Probability of success, p
266 static_cast<RealType>(1), // Probability of result (CDF), P
267 static_cast<RealType>(1.0000000000000364e-102), // Q = 1 - P
268 tolerance);
269
270 test_spot( // > formatC(pgeom(1,0.99, TRUE), digits=17) [1] "0.99990000000000001"
271 // > formatC(pgeom(1,0.99, FALSE), digits=17) [1] "0.00010000000000000009"
272 static_cast<RealType>(1), // Number of failures, k
273 static_cast<RealType>(0.99), // Probability of success, p
274 static_cast<RealType>(0.9999), // Probability of result (CDF), P
275 static_cast<RealType>(0.0001), // Q = 1 - P
276 tolerance);
277
278 if(std::numeric_limits<RealType>::is_specialized)
279 { // An extreme value test that is more accurate than using negative binomial.
280 // Since geometric only uses exp and log functions.
281 test_spot( // > formatC(pgeom(10000, 0.001, TRUE), digits=17) [1] "0.99995487182736897"
282 // > formatC(pgeom(10000,0.001, FALSE), digits=17) [1] "4.5128172631071587e-05"
283 static_cast<RealType>(10000L), // Number of failures, k
284 static_cast<RealType>(0.001L), // Probability of success, p
285 static_cast<RealType>(0.99995487182736897L), // Probability of result (CDF), P
286 static_cast<RealType>(4.5128172631071587e-05L), // Q = 1 - P
287 tolerance); //
288 } // numeric_limit is specialized
289 // End of single spot tests using RealType
290
291 // Tests on PDF:
292
293 BOOST_CHECK_CLOSE_FRACTION( //> formatC(dgeom(0,0.5), digits=17)[1] " 0.5"
294 pdf(geometric_distribution<RealType>(static_cast<RealType>(0.5)),
295 static_cast<RealType>(0.0) ), // Number of failures, k is very small but not integral,
296 static_cast<RealType>(0.5), // nearly success probability.
297 tolerance);
298
299 BOOST_CHECK_CLOSE_FRACTION( //> formatC(dgeom(0,0.5), digits=17)[1] " 0.5"
300 // R treates geom as a discrete distribution.
301 // > formatC(dgeom(1.999999,0.5, FALSE), digits=17) [1] " 0"
302 // Warning message:
303 // In dgeom(1.999999, 0.5, FALSE) : non-integer x = 1.999999
304 pdf(geometric_distribution<RealType>(static_cast<RealType>(0.5)),
305 static_cast<RealType>(0.0001L) ), // Number of failures, k is very small but not integral,
306 static_cast<RealType>(0.4999653438420768L), // nearly success probability.
307 tolerance);
308
309 BOOST_CHECK_CLOSE_FRACTION( // > formatC(pgeom(0.0001,0.5, TRUE), digits=17)[1] " 0.5"
310 // > formatC(pgeom(0.0001,0.5, FALSE), digits=17) [1] " 0.5"
311 // R treates geom as a discrete distribution.
312 pdf(geometric_distribution<RealType>(static_cast<RealType>(0.5)),
313 static_cast<RealType>(0.0001L) ), // Number of failures, k is very small but not integral,
314 static_cast<RealType>(0.4999653438420768L), // nearly success probability.
315 tolerance);
316
317 BOOST_CHECK_CLOSE_FRACTION( // formatC(dgeom(1,0.01), digits=17)[1] "0.0099000000000000008"
318 pdf(geometric_distribution<RealType>(static_cast<RealType>(0.01L)),
319 static_cast<RealType>(1) ), // Number of failures, k
320 static_cast<RealType>(0.0099000000000000008), //
321 tolerance);
322
323 BOOST_CHECK_CLOSE_FRACTION( //> formatC(dgeom(1,0.99), digits=17)[1] "0.0099000000000000043"
324 pdf(geometric_distribution<RealType>(static_cast<RealType>(0.99L)),
325 static_cast<RealType>(1) ), // Number of failures, k
326 static_cast<RealType>(0.00990000000000000043L), //
327 tolerance);
328
329 BOOST_CHECK_CLOSE_FRACTION( //> > formatC(dgeom(0,0.99), digits=17)[1] "0.98999999999999999"
330 pdf(geometric_distribution<RealType>(static_cast<RealType>(0.99L)),
331 static_cast<RealType>(0) ), // Number of failures, k
332 static_cast<RealType>(0.98999999999999999L), //
333 tolerance);
334
335 // p near unity.
336 BOOST_CHECK_CLOSE_FRACTION( // > formatC(dgeom(100,0.99), digits=17)[1] "9.9000000000003448e-201"
337 pdf(geometric_distribution<RealType>(static_cast<RealType>(0.99L)),
338 static_cast<RealType>(100) ), // Number of failures, k
339 static_cast<RealType>(9.9000000000003448e-201L), //
340 100 * tolerance); // Note difference
341
342 // p nearer unity.
343 BOOST_CHECK_CLOSE_FRACTION( //
344 pdf(geometric_distribution<RealType>(static_cast<RealType>(0.9999)),
345 static_cast<RealType>(10) ), // Number of failures, k
346 // static_cast<double>(9.9989999999889024e-41), // Boost.Math
347 // static_cast<float>(1.00156406e-040)
348 static_cast<RealType>(9.999e-41), // exact from 100 digit calculator.
349 2e3 * tolerance); // Note bigger tolerance needed.
350
351 // Moshier Cephes 100 digits calculator says 9.999e-41
352 //0.9999*pow(1-0.9999,10)
353 // 9.9990000000000000000000000000000000000000000000000000000000000000000000E-41
354 // 9.998999999988988e-041
355 // > formatC(dgeom(10, 0.9999), digits=17) [1] "9.9989999999889024e-41"
356 // p * pow(q, k) 9.9989999999889880e-041
357 // exp(p * k * log1p(-p)) 9.9989999999889024e-041
358
359
360
361 // 0.9999999999 * pow(1-0.9999999999,10)= 9.9999999990E-101
362 // > formatC(dgeom(10,0.9999999999), digits=17) [1] "1.0000008273040127e-100"
363 BOOST_CHECK_CLOSE_FRACTION( //
364 pdf(geometric_distribution<RealType>(static_cast<RealType>(0.9999999999L)),
365 static_cast<RealType>(10) ), //
366 static_cast<RealType>(9.9999999990E-101L), // 1.0000008273040179e-100
367 1e9 * tolerance); // Note big tolerance needed.
368 // 1.0000008273040179e-100 Boost.Math
369 // 1.0000008273040127e-100 R
370 // 0.9999999990000004e-100 100 digit calculator 'exact'
371
372 BOOST_CHECK_CLOSE_FRACTION( //
373 pdf(geometric_distribution<RealType>(static_cast<RealType>(0.00000000001L)),
374 static_cast<RealType>(10) ), //
375 static_cast<RealType>(9.999999999e-12L), // get 9.9999999989999994e-012
376 1 * tolerance); // Note small tolerance needed.
377
378
379 BOOST_CHECK_CLOSE_FRACTION( //
380 pdf(geometric_distribution<RealType>(static_cast<RealType>(0.00000000001L)),
381 static_cast<RealType>(1000) ), //
382 static_cast<RealType>(9.9999999e-12L), // get 9.9999998999999913e-012
383 tolerance); // Note small tolerance needed.
384
385
386 ///////////////////////////////////////////////////
387 BOOST_CHECK_CLOSE_FRACTION( //
388 // > formatC(dgeom(0.0001,0.5, FALSE), digits=17) [1] " 0.5"
389 // R treates geom as a discrete distribution.
390 // But Boost.Math is continuous, so if you want R behaviour,
391 // make number of failures, k into an integer with the floor function.
392 pdf(geometric_distribution<RealType>(static_cast<RealType>(0.5)),
393 static_cast<RealType>(floor(0.0001L)) ), // Number of failures, k is very small but MADE integral,
394 static_cast<RealType>(0.5), // nearly success probability.
395 tolerance);
396
397 // R switches over at about 1e7 from k = 0, returning 0.5, to k = 1, returning 0.25.
398 // Boost.Math does not do this, even for 0.9999999999999999
399 // > formatC(pgeom(0.999999,0.5, FALSE), digits=17) [1] " 0.5"
400 // > formatC(pgeom(0.9999999,0.5, FALSE), digits=17) [1] " 0.25"
401
402 BOOST_CHECK_CLOSE_FRACTION( // > formatC(pgeom(0.0001,0.5, TRUE), digits=17)[1] " 0.5"
403 // > formatC(pgeom(0.0001,0.5, FALSE), digits=17) [1] " 0.5"
404 // R treates geom as a discrete distribution.
405 // But Boost.Math is continuous, so if you want R behaviour,
406 // make number of failures, k into an integer with the floor function.
407 pdf(geometric_distribution<RealType>(static_cast<RealType>(0.5)),
408 static_cast<RealType>(floor(0.9999999999999999L)) ), // Number of failures, k is very small but MADE integral,
409 static_cast<RealType>(0.5), // nearly success probability.
410 tolerance);
411
412 BOOST_CHECK_CLOSE_FRACTION( // > formatC(pgeom(0.0001,0.5, TRUE), digits=17)[1] " 0.5"
413 // > formatC(pgeom(0.0001,0.5, FALSE), digits=17) [1] " 0.5"
414 // R treates geom as a discrete distribution.
415 // But Boost.Math is continuous, so if you want R behaviour,
416 // make number of failures, k into an integer with the floor function.
417 pdf(geometric_distribution<RealType>(static_cast<RealType>(0.5)),
418 static_cast<RealType>(floor(1. - tolerance)) ),
419 // Number of failures, k is very small but MADE integral,
420 // Need to use tolerance here,
421 // as epsilon is ill-defined for Real concept:
422 // numeric_limits<RealType>::epsilon() 0
423 static_cast<RealType>(0.5), // nearly success probability.
424 tolerance * 10);
425
426 BOOST_CHECK_CLOSE_FRACTION(
427 pdf(geometric_distribution<RealType>(static_cast<RealType>(0.0001L)),
428 static_cast<RealType>(2)), // k = 2.
429 static_cast<RealType>(9.99800010e-5L), // 'exact '
430 tolerance);
431
432 //> formatC(dgeom(2, 0.9999), digits=17) [1] "9.9989999999977806e-09"
433 BOOST_CHECK_CLOSE_FRACTION(
434 pdf(geometric_distribution<RealType>(static_cast<RealType>(0.9999L)),
435 static_cast<RealType>(2)), // k = 0
436 static_cast<RealType>(9.999e-9L), // 'exact'
437 1000*tolerance);
438
439 BOOST_CHECK_CLOSE_FRACTION(
440 pdf(geometric_distribution<RealType>(static_cast<RealType>(0.9999L)),
441 static_cast<RealType>(3)), // k = 3
442 static_cast<RealType>(9.999e-13L), // get
443 1000*tolerance);
444
445 BOOST_CHECK_CLOSE_FRACTION(
446 pdf(geometric_distribution<RealType>(static_cast<RealType>(0.9999L)),
447 static_cast<RealType>(5)), // k = 5
448 static_cast<RealType>(9.999e-21L), // 9.9989999999944947e-021
449 1000*tolerance);
450
451
452 BOOST_CHECK_CLOSE_FRACTION(
453 pdf(geometric_distribution<RealType>( static_cast<RealType>(0.0001L)),
454 static_cast<RealType>(3)), // k = 0.
455 static_cast<RealType>(9.99700029999e-5L), //
456 tolerance);
457 // Tests on cdf:
458 // MathCAD pgeom k, r, p) == failures, successes, probability.
459
460 BOOST_CHECK_CLOSE_FRACTION(cdf(
461 geometric_distribution<RealType>(static_cast<RealType>(0.5)), // prob 0.5
462 static_cast<RealType>(0) ), // k = 0
463 static_cast<RealType>(0.5), // probability =p
464 tolerance);
465
466 BOOST_CHECK_CLOSE_FRACTION(cdf(complement(
467 geometric_distribution<RealType>(static_cast<RealType>(0.5)), //
468 static_cast<RealType>(0) )), // k = 0
469 static_cast<RealType>(0.5), // probability =
470 tolerance);
471
472 BOOST_CHECK_CLOSE_FRACTION(cdf(
473 geometric_distribution<RealType>(static_cast<RealType>(0.25)), // prob 0.5
474 static_cast<RealType>(1) ), // k = 0
475 static_cast<RealType>(0.4375L), // probability =p
476 tolerance);
477
478 BOOST_CHECK_CLOSE_FRACTION(cdf(complement(
479 geometric_distribution<RealType>(static_cast<RealType>(0.25)), //
480 static_cast<RealType>(1) )), // k = 0
481 static_cast<RealType>(1-0.4375L), // probability =
482 tolerance);
483
484 BOOST_CHECK_CLOSE_FRACTION(cdf(complement(
485 geometric_distribution<RealType>(static_cast<RealType>(0.5)), //
486 static_cast<RealType>(1) )), // k = 0
487 static_cast<RealType>(0.25), // probability = exact 0.25
488 tolerance);
489
490 BOOST_CHECK_CLOSE_FRACTION( //
491 cdf(geometric_distribution<RealType>(static_cast<RealType>(0.5)),
492 static_cast<RealType>(4)), // k =4.
493 static_cast<RealType>(0.96875L), // exact
494 tolerance);
495
496
497 // Tests of other functions, mean and other moments ...
498
499 geometric_distribution<RealType> dist(static_cast<RealType>(0.25));
500 // mean:
501 BOOST_CHECK_CLOSE_FRACTION(
502 mean(dist), static_cast<RealType>((1 - 0.25) /0.25), tol5eps);
503 BOOST_CHECK_CLOSE_FRACTION(
504 mode(dist), static_cast<RealType>(0), tol1eps);
505 // variance:
506 BOOST_CHECK_CLOSE_FRACTION(
507 variance(dist), static_cast<RealType>((1 - 0.25) / (0.25 * 0.25)), tol5eps);
508
509 // std deviation:
510 // sqrt(0.75/0.125)
511
512 BOOST_CHECK_CLOSE_FRACTION(
513 standard_deviation(dist), //
514 static_cast<RealType>(sqrt((1.0L - 0.25L) / (0.25L * 0.25L))), // using 100 digit calc
515 tol5eps);
516
517 BOOST_CHECK_CLOSE_FRACTION(
518 skewness(dist), //
519 static_cast<RealType>((2-0.25L) /sqrt(0.75L)),
520 // using calculator
521 tol5eps);
522 BOOST_CHECK_CLOSE_FRACTION(
523 kurtosis_excess(dist), //
524 static_cast<RealType>(6 + 0.0625L/0.75L), //
525 tol5eps);
526 // 6.083333333333333 6.166666666666667
527 BOOST_CHECK_CLOSE_FRACTION(
528 kurtosis(dist), // true
529 static_cast<RealType>(9 + 0.0625L/0.75L), //
530 tol5eps);
531 // hazard:
532 RealType x = static_cast<RealType>(0.125);
533 BOOST_CHECK_CLOSE_FRACTION(
534 hazard(dist, x)
535 , pdf(dist, x) / cdf(complement(dist, x)), tol5eps);
536 // cumulative hazard:
537 BOOST_CHECK_CLOSE_FRACTION(
538 chf(dist, x), -log(cdf(complement(dist, x))), tol5eps);
539 // coefficient_of_variation:
540 BOOST_CHECK_CLOSE_FRACTION(
541 coefficient_of_variation(dist)
542 , standard_deviation(dist) / mean(dist), tol5eps);
543
544 // Special cases for PDF:
545 BOOST_CHECK_EQUAL(
546 pdf(
547 geometric_distribution<RealType>(static_cast<RealType>(0)), //
548 static_cast<RealType>(0)),
549 static_cast<RealType>(0) );
550
551 BOOST_CHECK_EQUAL(
552 pdf(
553 geometric_distribution<RealType>(static_cast<RealType>(0)),
554 static_cast<RealType>(0.0001)),
555 static_cast<RealType>(0) );
556
557 BOOST_CHECK_EQUAL(
558 pdf(
559 geometric_distribution<RealType>(static_cast<RealType>(1)),
560 static_cast<RealType>(0.001)),
561 static_cast<RealType>(0) );
562
563 BOOST_CHECK_EQUAL(
564 pdf(
565 geometric_distribution<RealType>(static_cast<RealType>(1)),
566 static_cast<RealType>(8)),
567 static_cast<RealType>(0) );
568
569 BOOST_CHECK_SMALL(
570 pdf(
571 geometric_distribution<RealType>(static_cast<RealType>(0.25)),
572 static_cast<RealType>(0))-
573 static_cast<RealType>(0.25),
574 2 * boost::math::tools::epsilon<RealType>() ); // Expect exact, but not quite.
575 // numeric_limits<RealType>::epsilon()); // Not suitable for real concept!
576
577 // Quantile boundary cases checks:
578 BOOST_CHECK_EQUAL(
579 quantile( // zero P < cdf(0) so should be exactly zero.
580 geometric_distribution<RealType>(static_cast<RealType>(0.25)),
581 static_cast<RealType>(0)),
582 static_cast<RealType>(0));
583
584 BOOST_CHECK_EQUAL(
585 quantile( // min P < cdf(0) so should be exactly zero.
586 geometric_distribution<RealType>(static_cast<RealType>(0.25)),
587 static_cast<RealType>(boost::math::tools::min_value<RealType>())),
588 static_cast<RealType>(0));
589
590 BOOST_CHECK_CLOSE_FRACTION(
591 quantile( // Small P < cdf(0) so should be near zero.
592 geometric_distribution<RealType>(static_cast<RealType>(0.25)),
593 static_cast<RealType>(boost::math::tools::epsilon<RealType>())), //
594 static_cast<RealType>(0),
595 tol5eps);
596
597 BOOST_CHECK_CLOSE_FRACTION(
598 quantile( // Small P < cdf(0) so should be exactly zero.
599 geometric_distribution<RealType>(static_cast<RealType>(0.25)),
600 static_cast<RealType>(0.0001)),
601 static_cast<RealType>(0),
602 tolerance);
603
604 //BOOST_CHECK( // Fails with overflow for real_concept
605 //quantile( // Small P near 1 so k failures should be big.
606 //geometric_distribution<RealType>(static_cast<RealType>(8), static_cast<RealType>(0.25)),
607 //static_cast<RealType>(1 - boost::math::tools::epsilon<RealType>())) <=
608 //static_cast<RealType>(189.56999032670058) // 106.462769 for float
609 //);
610
611 if(std::numeric_limits<RealType>::has_infinity)
612 { // BOOST_CHECK tests for infinity using std::numeric_limits<>::infinity()
613 // Note that infinity is not implemented for real_concept, so these tests
614 // are only done for types, like built-in float, double.. that have infinity.
615 // Note that these assume that BOOST_MATH_OVERFLOW_ERROR_POLICY is NOT throw_on_error.
616 // #define BOOST_MATH_THROW_ON_OVERFLOW_POLICY == throw_on_error would throw here.
617 // #define BOOST_MAT_DOMAIN_ERROR_POLICY IS defined throw_on_error,
618 // so the throw path of error handling is tested below with BOOST_CHECK_THROW tests.
619
620 BOOST_CHECK(
621 quantile( // At P == 1 so k failures should be infinite.
622 geometric_distribution<RealType>(static_cast<RealType>(0.25)),
623 static_cast<RealType>(1)) ==
624 //static_cast<RealType>(boost::math::tools::infinity<RealType>())
625 static_cast<RealType>(std::numeric_limits<RealType>::infinity()) );
626
627 BOOST_CHECK_EQUAL(
628 quantile( // At 1 == P so should be infinite.
629 geometric_distribution<RealType>( static_cast<RealType>(0.25)),
630 static_cast<RealType>(1)), //
631 std::numeric_limits<RealType>::infinity() );
632
633 BOOST_CHECK_EQUAL(
634 quantile(complement( // Q zero 1 so P == 1 < cdf(0) so should be exactly infinity.
635 geometric_distribution<RealType>(static_cast<RealType>(0.25)),
636 static_cast<RealType>(0))),
637 std::numeric_limits<RealType>::infinity() );
638 } // test for infinity using std::numeric_limits<>::infinity()
639 else
640 { // real_concept case, so check it throws rather than returning infinity.
641 BOOST_CHECK_EQUAL(
642 quantile( // At P == 1 so k failures should be infinite.
643 geometric_distribution<RealType>(static_cast<RealType>(0.25)),
644 static_cast<RealType>(1)),
645 boost::math::tools::max_value<RealType>() );
646
647 BOOST_CHECK_EQUAL(
648 quantile(complement( // Q zero 1 so P == 1 < cdf(0) so should be exactly infinity.
649 geometric_distribution<RealType>(static_cast<RealType>(0.25)),
650 static_cast<RealType>(0))),
651 boost::math::tools::max_value<RealType>());
652 } // has infinity
653
654 BOOST_CHECK( // Should work for built-in and real_concept.
655 quantile(complement( // Q near to 1 so P nearly 1, so should be large > 300.
656 geometric_distribution<RealType>(static_cast<RealType>(0.25)),
657 static_cast<RealType>(boost::math::tools::min_value<RealType>())))
658 >= static_cast<RealType>(300) );
659
660 BOOST_CHECK_EQUAL(
661 quantile( // P == 0 < cdf(0) so should be zero.
662 geometric_distribution<RealType>(static_cast<RealType>(0.25)),
663 static_cast<RealType>(0)),
664 static_cast<RealType>(0));
665
666 // Quantile Complement boundary cases:
667
668 BOOST_CHECK_EQUAL(
669 quantile(complement( // Q = 1 so P = 0 < cdf(0) so should be exactly zero.
670 geometric_distribution<RealType>( static_cast<RealType>(0.25)),
671 static_cast<RealType>(1))),
672 static_cast<RealType>(0)
673 );
674
675 BOOST_CHECK_EQUAL(
676 quantile(complement( // Q very near 1 so P == epsilon < cdf(0) so should be exactly zero.
677 geometric_distribution<RealType>(static_cast<RealType>(0.25)),
678 static_cast<RealType>(1 - boost::math::tools::epsilon<RealType>()))),
679 static_cast<RealType>(0)
680 );
681
682 // Check that duff arguments throw domain_error:
683
684 BOOST_CHECK_THROW(
685 pdf( // Negative success_fraction!
686 geometric_distribution<RealType>(static_cast<RealType>(-0.25)),
687 static_cast<RealType>(0)), std::domain_error);
688 BOOST_CHECK_THROW(
689 pdf( // Success_fraction > 1!
690 geometric_distribution<RealType>(static_cast<RealType>(1.25)),
691 static_cast<RealType>(0)),
692 std::domain_error);
693 BOOST_CHECK_THROW(
694 pdf( // Negative k argument !
695 geometric_distribution<RealType>(static_cast<RealType>(0.25)),
696 static_cast<RealType>(-1)),
697 std::domain_error);
698 //BOOST_CHECK_THROW(
699 //pdf( // check limit on k (failures)
700 //geometric_distribution<RealType>(static_cast<RealType>(0.25)),
701 //std::numeric_limits<RealType>infinity()),
702 //std::domain_error);
703 BOOST_CHECK_THROW(
704 cdf( // Negative k argument !
705 geometric_distribution<RealType>(static_cast<RealType>(0.25)),
706 static_cast<RealType>(-1)),
707 std::domain_error);
708 BOOST_CHECK_THROW(
709 cdf( // Negative success_fraction!
710 geometric_distribution<RealType>(static_cast<RealType>(-0.25)),
711 static_cast<RealType>(0)), std::domain_error);
712 BOOST_CHECK_THROW(
713 cdf( // Success_fraction > 1!
714 geometric_distribution<RealType>(static_cast<RealType>(1.25)),
715 static_cast<RealType>(0)), std::domain_error);
716 BOOST_CHECK_THROW(
717 quantile( // Negative success_fraction!
718 geometric_distribution<RealType>(static_cast<RealType>(-0.25)),
719 static_cast<RealType>(0)), std::domain_error);
720 BOOST_CHECK_THROW(
721 quantile( // Success_fraction > 1!
722 geometric_distribution<RealType>(static_cast<RealType>(1.25)),
723 static_cast<RealType>(0)), std::domain_error);
724 check_out_of_range<geometric_distribution<RealType> >(0.5);
725 // End of check throwing 'duff' out-of-domain values.
726
727 { // Compare geometric and negative binomial functions.
728 using boost::math::negative_binomial_distribution;
729 using boost::math::geometric_distribution;
730
731 RealType k = static_cast<RealType>(2.L);
732 RealType alpha = static_cast<RealType>(0.05L);
733 RealType p = static_cast<RealType>(0.5L);
734
735 BOOST_CHECK_CLOSE_FRACTION( // Successes parameter in negative binomial is 1 for geometric.
736 geometric_distribution<RealType>::find_lower_bound_on_p(k, alpha),
737 negative_binomial_distribution<RealType>::find_lower_bound_on_p(k, static_cast<RealType>(1), alpha),
738 tolerance);
739 BOOST_CHECK_CLOSE_FRACTION( // Successes parameter in negative binomial is 1 for geometric.
740 geometric_distribution<RealType>::find_upper_bound_on_p(k, alpha),
741 negative_binomial_distribution<RealType>::find_upper_bound_on_p(k, static_cast<RealType>(1), alpha),
742 tolerance);
743 BOOST_CHECK_CLOSE_FRACTION( // Should be identical - successes parameter is not used.
744 geometric_distribution<RealType>::find_maximum_number_of_trials(k, p, alpha),
745 negative_binomial_distribution<RealType>::find_maximum_number_of_trials(k, p, alpha),
746 tolerance);
747 }
748 //geometric::find_upper_bound_on_p(k, alpha);
749 return;
750 } // template <class RealType> void test_spots(RealType) // Any floating-point type RealType.
751
BOOST_AUTO_TEST_CASE(test_main)752 BOOST_AUTO_TEST_CASE( test_main )
753 {
754 // Check that can generate geometric distribution using the two convenience methods:
755 using namespace boost::math;
756 geometric g05d(0.5); // Using typedef - default type is double.
757 geometric_distribution<> g05dd(0.5); // Using default RealType double.
758
759 // Basic sanity-check spot values.
760
761 // Test some simple double only examples.
762 geometric_distribution<double> mydist(0.25);
763 // success fraction == 0.25 == 25% or 1 in 4 successes.
764 // Note: double values (matching the distribution definition) avoid the need for any casting.
765
766 // Check accessor functions return exact values for double at least.
767 BOOST_CHECK_EQUAL(mydist.success_fraction(), static_cast<double>(1./4.));
768
769 //cout << numeric_limits<RealType>::epsilon() << endl;
770
771 // (Parameter value, arbitrarily zero, only communicates the floating point type).
772 #ifdef TEST_FLOAT
773 test_spots(0.0F); // Test float.
774 #endif
775 #ifdef TEST_DOUBLE
776 test_spots(0.0); // Test double.
777 #endif
778 #ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
779 #ifdef TEST_LDOUBLE
780 test_spots(0.0L); // Test long double.
781 #endif
782 #if !BOOST_WORKAROUND(__BORLANDC__, BOOST_TESTED_AT(0x582))
783 #ifdef TEST_REAL_CONCEPT
784 test_spots(boost::math::concepts::real_concept(0.)); // Test real concept.
785 #endif
786 #endif
787 #else
788 std::cout << "<note>The long double tests have been disabled on this platform "
789 "either because the long double overloads of the usual math functions are "
790 "not available at all, or because they are too inaccurate for these tests "
791 "to pass.</note>" << std::cout;
792 #endif
793
794
795 } // BOOST_AUTO_TEST_CASE( test_main )
796
797 /*
798
799
800
801 */
802