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