]>
Commit | Line | Data |
---|---|---|
7c673cae FG |
1 | // test_binomial.cpp |
2 | ||
3 | // Copyright John Maddock 2006. | |
4 | // Copyright Paul A. Bristow 2007. | |
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 | // Basic sanity test for Binomial Cumulative Distribution Function. | |
12 | ||
13 | #define BOOST_MATH_DISCRETE_QUANTILE_POLICY real | |
14 | ||
15 | #if !defined(TEST_FLOAT) && !defined(TEST_DOUBLE) && !defined(TEST_LDOUBLE) && !defined(TEST_REAL_CONCEPT) | |
16 | # define TEST_FLOAT | |
17 | # define TEST_DOUBLE | |
18 | # define TEST_LDOUBLE | |
19 | # define TEST_REAL_CONCEPT | |
20 | #endif | |
21 | ||
22 | #ifdef _MSC_VER | |
23 | # pragma warning(disable: 4127) // conditional expression is constant. | |
24 | # pragma warning(disable: 4100) // unreferenced formal parameter. | |
25 | // Seems an entirely spurious warning - formal parameter T IS used - get error if /* T */ | |
26 | //# pragma warning(disable: 4535) // calling _set_se_translator() requires /EHa (in Boost.test) | |
27 | // Enable C++ Exceptions Yes With SEH Exceptions (/EHa) prevents warning 4535. | |
28 | #endif | |
29 | ||
30 | #include <boost/math/tools/test.hpp> | |
31 | #include <boost/math/concepts/real_concept.hpp> // for real_concept | |
32 | using ::boost::math::concepts::real_concept; | |
33 | ||
34 | #include <boost/math/distributions/binomial.hpp> // for binomial_distribution | |
35 | using boost::math::binomial_distribution; | |
36 | ||
37 | #define BOOST_TEST_MAIN | |
38 | #include <boost/test/unit_test.hpp> // for test_main | |
39 | #include <boost/test/floating_point_comparison.hpp> // for BOOST_CHECK_CLOSE | |
40 | #include "table_type.hpp" | |
41 | ||
42 | #include "test_out_of_range.hpp" | |
43 | ||
44 | #include <iostream> | |
45 | using std::cout; | |
46 | using std::endl; | |
47 | #include <limits> | |
48 | using std::numeric_limits; | |
49 | ||
50 | template <class RealType> | |
51 | void test_spot( | |
52 | RealType N, // Number of trials | |
53 | RealType k, // Number of successes | |
54 | RealType p, // Probability of success | |
55 | RealType P, // CDF | |
56 | RealType Q, // Complement of CDF | |
57 | RealType tol) // Test tolerance | |
58 | { | |
59 | boost::math::binomial_distribution<RealType> bn(N, p); | |
60 | BOOST_CHECK_CLOSE( | |
61 | cdf(bn, k), P, tol); | |
62 | if((P < 0.99) && (Q < 0.99)) | |
63 | { | |
64 | // | |
65 | // We can only check this if P is not too close to 1, | |
66 | // so that we can guarantee Q is free of error: | |
67 | // | |
68 | BOOST_CHECK_CLOSE( | |
69 | cdf(complement(bn, k)), Q, tol); | |
70 | if(k != 0) | |
71 | { | |
72 | BOOST_CHECK_CLOSE( | |
73 | quantile(bn, 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) && (boost::is_floating_point<RealType>::value)) | |
79 | { | |
80 | // Limit where this is checked: if exponent range is very large we may | |
81 | // run out of iterations in our root finding algorithm. | |
82 | BOOST_CHECK(quantile(bn, P) < boost::math::tools::epsilon<RealType>() * 10); | |
83 | } | |
84 | } | |
85 | if(k != 0) | |
86 | { | |
87 | BOOST_CHECK_CLOSE( | |
88 | quantile(complement(bn, Q)), k, tol); | |
89 | } | |
90 | else | |
91 | { | |
92 | // Just check quantile is very small: | |
93 | if((std::numeric_limits<RealType>::max_exponent <= std::numeric_limits<double>::max_exponent) && (boost::is_floating_point<RealType>::value)) | |
94 | { | |
95 | // Limit where this is checked: if exponent range is very large we may | |
96 | // run out of iterations in our root finding algorithm. | |
97 | BOOST_CHECK(quantile(complement(bn, Q)) < boost::math::tools::epsilon<RealType>() * 10); | |
98 | } | |
99 | } | |
100 | if(k > 0) | |
101 | { | |
102 | // estimate success ratio: | |
103 | // Note lower bound uses a different formual internally | |
104 | // from upper bound, have to adjust things to prevent | |
105 | // fencepost errors: | |
106 | BOOST_CHECK_CLOSE( | |
107 | binomial_distribution<RealType>::find_lower_bound_on_p( | |
108 | N, k+1, Q), | |
109 | p, tol); | |
110 | BOOST_CHECK_CLOSE( | |
111 | binomial_distribution<RealType>::find_upper_bound_on_p( | |
112 | N, k, P), | |
113 | p, tol); | |
114 | ||
115 | if(Q < P) | |
116 | { | |
117 | // Default method (Clopper Pearson) | |
118 | BOOST_CHECK( | |
119 | binomial_distribution<RealType>::find_lower_bound_on_p( | |
120 | N, k, Q) | |
121 | <= | |
122 | binomial_distribution<RealType>::find_upper_bound_on_p( | |
123 | N, k, Q) | |
124 | ); | |
125 | BOOST_CHECK(( | |
126 | binomial_distribution<RealType>::find_lower_bound_on_p( | |
127 | N, k, Q) | |
128 | <= k/N) && (k/N <= | |
129 | binomial_distribution<RealType>::find_upper_bound_on_p( | |
130 | N, k, Q)) | |
131 | ); | |
132 | // Bayes Method (Jeffreys Prior) | |
133 | BOOST_CHECK( | |
134 | binomial_distribution<RealType>::find_lower_bound_on_p( | |
135 | N, k, Q, binomial_distribution<RealType>::jeffreys_prior_interval) | |
136 | <= | |
137 | binomial_distribution<RealType>::find_upper_bound_on_p( | |
138 | N, k, Q, binomial_distribution<RealType>::jeffreys_prior_interval) | |
139 | ); | |
140 | BOOST_CHECK(( | |
141 | binomial_distribution<RealType>::find_lower_bound_on_p( | |
142 | N, k, Q, binomial_distribution<RealType>::jeffreys_prior_interval) | |
143 | <= k/N) && (k/N <= | |
144 | binomial_distribution<RealType>::find_upper_bound_on_p( | |
145 | N, k, Q, binomial_distribution<RealType>::jeffreys_prior_interval)) | |
146 | ); | |
147 | } | |
148 | else | |
149 | { | |
150 | // Default method (Clopper Pearson) | |
151 | BOOST_CHECK( | |
152 | binomial_distribution<RealType>::find_lower_bound_on_p( | |
153 | N, k, P) | |
154 | <= | |
155 | binomial_distribution<RealType>::find_upper_bound_on_p( | |
156 | N, k, P) | |
157 | ); | |
158 | BOOST_CHECK( | |
159 | (binomial_distribution<RealType>::find_lower_bound_on_p( | |
160 | N, k, P) | |
161 | <= k / N) && (k/N <= | |
162 | binomial_distribution<RealType>::find_upper_bound_on_p( | |
163 | N, k, P)) | |
164 | ); | |
165 | // Bayes Method (Jeffreys Prior) | |
166 | BOOST_CHECK( | |
167 | binomial_distribution<RealType>::find_lower_bound_on_p( | |
168 | N, k, P, binomial_distribution<RealType>::jeffreys_prior_interval) | |
169 | <= | |
170 | binomial_distribution<RealType>::find_upper_bound_on_p( | |
171 | N, k, P, binomial_distribution<RealType>::jeffreys_prior_interval) | |
172 | ); | |
173 | BOOST_CHECK( | |
174 | (binomial_distribution<RealType>::find_lower_bound_on_p( | |
175 | N, k, P, binomial_distribution<RealType>::jeffreys_prior_interval) | |
176 | <= k / N) && (k/N <= | |
177 | binomial_distribution<RealType>::find_upper_bound_on_p( | |
178 | N, k, P, binomial_distribution<RealType>::jeffreys_prior_interval)) | |
179 | ); | |
180 | } | |
181 | } | |
182 | // | |
183 | // estimate sample size: | |
184 | // | |
185 | BOOST_CHECK_CLOSE( | |
186 | binomial_distribution<RealType>::find_minimum_number_of_trials( | |
187 | k, p, P), | |
188 | N, tol); | |
189 | BOOST_CHECK_CLOSE( | |
190 | binomial_distribution<RealType>::find_maximum_number_of_trials( | |
191 | k, p, Q), | |
192 | N, tol); | |
193 | } | |
194 | ||
195 | // Double check consistency of CDF and PDF by computing | |
196 | // the finite sum: | |
197 | RealType sum = 0; | |
198 | for(unsigned i = 0; i <= k; ++i) | |
199 | sum += pdf(bn, RealType(i)); | |
200 | BOOST_CHECK_CLOSE( | |
201 | sum, P, tol); | |
202 | // And complement as well: | |
203 | sum = 0; | |
204 | for(RealType i = N; i > k; i -= 1) | |
205 | sum += pdf(bn, i); | |
206 | if(P < 0.99) | |
207 | { | |
208 | BOOST_CHECK_CLOSE( | |
209 | sum, Q, tol); | |
210 | } | |
211 | else | |
212 | { | |
213 | // Not enough information content in P for Q to be meaningful | |
214 | RealType tol = (std::max)(2 * Q, boost::math::tools::epsilon<RealType>()); | |
215 | BOOST_CHECK(sum < tol); | |
216 | } | |
217 | } | |
218 | ||
219 | template <class RealType> // Any floating-point type RealType. | |
220 | void test_spots(RealType T) | |
221 | { | |
222 | // Basic sanity checks, test data is to double precision only | |
223 | // so set tolerance to 100eps expressed as a persent, or | |
224 | // 100eps of type double expressed as a persent, whichever | |
225 | // is the larger. | |
226 | ||
227 | RealType tolerance = (std::max) | |
228 | (boost::math::tools::epsilon<RealType>(), | |
229 | static_cast<RealType>(std::numeric_limits<double>::epsilon())); | |
230 | tolerance *= 100 * 1000; | |
231 | RealType tol2 = boost::math::tools::epsilon<RealType>() * 5 * 100; // 5 eps as a persent | |
232 | ||
233 | cout << "Tolerance for type " << typeid(T).name() << " is " << tolerance << " %" << endl; | |
234 | ||
235 | ||
236 | // Sources of spot test values: | |
237 | ||
238 | // MathCAD defines pbinom(k, n, p) | |
239 | // returns pr(X ,=k) when random variable X has the binomial distribution with parameters n and p. | |
240 | // 0 <= k ,= n | |
241 | // 0 <= p <= 1 | |
242 | // P = pbinom(30, 500, 0.05) = 0.869147702104609 | |
243 | ||
244 | using boost::math::binomial_distribution; | |
245 | using ::boost::math::cdf; | |
246 | using ::boost::math::pdf; | |
247 | ||
248 | #if !defined(TEST_ROUNDING) || (TEST_ROUNDING == 0) | |
249 | // Test binomial using cdf spot values from MathCAD. | |
250 | // These test quantiles and complements as well. | |
251 | test_spot( | |
252 | static_cast<RealType>(500), // Sample size, N | |
253 | static_cast<RealType>(30), // Number of successes, k | |
254 | static_cast<RealType>(0.05), // Probability of success, p | |
255 | static_cast<RealType>(0.869147702104609), // Probability of result (CDF), P | |
256 | static_cast<RealType>(1 - 0.869147702104609), // Q = 1 - P | |
257 | tolerance); | |
258 | ||
259 | test_spot( | |
260 | static_cast<RealType>(500), // Sample size, N | |
261 | static_cast<RealType>(250), // Number of successes, k | |
262 | static_cast<RealType>(0.05), // Probability of success, p | |
263 | static_cast<RealType>(1), // Probability of result (CDF), P | |
264 | static_cast<RealType>(0), // Q = 1 - P | |
265 | tolerance); | |
266 | ||
267 | test_spot( | |
268 | static_cast<RealType>(500), // Sample size, N | |
269 | static_cast<RealType>(470), // Number of successes, k | |
270 | static_cast<RealType>(0.95), // Probability of success, p | |
271 | static_cast<RealType>(0.176470742656766), // Probability of result (CDF), P | |
272 | static_cast<RealType>(1 - 0.176470742656766), // Q = 1 - P | |
273 | tolerance * 10); // Note higher tolerance on this test! | |
274 | ||
275 | test_spot( | |
276 | static_cast<RealType>(500), // Sample size, N | |
277 | static_cast<RealType>(400), // Number of successes, k | |
278 | static_cast<RealType>(0.05), // Probability of success, p | |
279 | static_cast<RealType>(1), // Probability of result (CDF), P | |
280 | static_cast<RealType>(0), // Q = 1 - P | |
281 | tolerance); | |
282 | ||
283 | test_spot( | |
284 | static_cast<RealType>(500), // Sample size, N | |
285 | static_cast<RealType>(400), // Number of successes, k | |
286 | static_cast<RealType>(0.9), // Probability of success, p | |
287 | static_cast<RealType>(1.80180425681923E-11), // Probability of result (CDF), P | |
288 | static_cast<RealType>(1 - 1.80180425681923E-11), // Q = 1 - P | |
289 | tolerance); | |
290 | ||
291 | test_spot( | |
292 | static_cast<RealType>(500), // Sample size, N | |
293 | static_cast<RealType>(5), // Number of successes, k | |
294 | static_cast<RealType>(0.05), // Probability of success, p | |
295 | static_cast<RealType>(9.181808267643E-7), // Probability of result (CDF), P | |
296 | static_cast<RealType>(1 - 9.181808267643E-7), // Q = 1 - P | |
297 | tolerance); | |
298 | ||
299 | test_spot( | |
300 | static_cast<RealType>(2), // Sample size, N | |
301 | static_cast<RealType>(1), // Number of successes, k | |
302 | static_cast<RealType>(0.5), // Probability of success, p | |
303 | static_cast<RealType>(0.75), // Probability of result (CDF), P | |
304 | static_cast<RealType>(0.25), // Q = 1 - P | |
305 | tolerance); | |
306 | ||
307 | test_spot( | |
308 | static_cast<RealType>(8), // Sample size, N | |
309 | static_cast<RealType>(3), // Number of successes, k | |
310 | static_cast<RealType>(0.25), // Probability of success, p | |
311 | static_cast<RealType>(0.8861846923828125), // Probability of result (CDF), P | |
312 | static_cast<RealType>(1 - 0.8861846923828125), // Q = 1 - P | |
313 | tolerance); | |
314 | ||
315 | test_spot( | |
316 | static_cast<RealType>(8), // Sample size, N | |
317 | static_cast<RealType>(0), // Number of successes, k | |
318 | static_cast<RealType>(0.25), // Probability of success, p | |
319 | static_cast<RealType>(0.1001129150390625), // Probability of result (CDF), P | |
320 | static_cast<RealType>(1 - 0.1001129150390625), // Q = 1 - P | |
321 | tolerance); | |
322 | ||
323 | test_spot( | |
324 | static_cast<RealType>(8), // Sample size, N | |
325 | static_cast<RealType>(1), // Number of successes, k | |
326 | static_cast<RealType>(0.25), // Probability of success, p | |
327 | static_cast<RealType>(0.36708068847656244), // Probability of result (CDF), P | |
328 | static_cast<RealType>(1 - 0.36708068847656244), // Q = 1 - P | |
329 | tolerance); | |
330 | ||
331 | test_spot( | |
332 | static_cast<RealType>(8), // Sample size, N | |
333 | static_cast<RealType>(4), // Number of successes, k | |
334 | static_cast<RealType>(0.25), // Probability of success, p | |
335 | static_cast<RealType>(0.9727020263671875), // Probability of result (CDF), P | |
336 | static_cast<RealType>(1 - 0.9727020263671875), // Q = 1 - P | |
337 | tolerance); | |
338 | ||
339 | test_spot( | |
340 | static_cast<RealType>(8), // Sample size, N | |
341 | static_cast<RealType>(7), // Number of successes, k | |
342 | static_cast<RealType>(0.25), // Probability of success, p | |
343 | static_cast<RealType>(0.9999847412109375), // Probability of result (CDF), P | |
344 | static_cast<RealType>(1 - 0.9999847412109375), // Q = 1 - P | |
345 | tolerance); | |
346 | ||
347 | // Tests on PDF follow: | |
348 | BOOST_CHECK_CLOSE( | |
349 | pdf(binomial_distribution<RealType>(static_cast<RealType>(20), static_cast<RealType>(0.75)), | |
350 | static_cast<RealType>(10)), // k. | |
351 | static_cast<RealType>(0.00992227527967770583927631378173), // 0.00992227527967770583927631378173 | |
352 | tolerance); | |
353 | ||
354 | BOOST_CHECK_CLOSE( | |
355 | pdf(binomial_distribution<RealType>(static_cast<RealType>(20), static_cast<RealType>(0.5)), | |
356 | static_cast<RealType>(10)), // k. | |
357 | static_cast<RealType>(0.17619705200195312500000000000000000000), // get k=10 0.049611376398388612 p = 0.25 | |
358 | tolerance); | |
359 | ||
360 | // Binomial pdf Test values from | |
361 | // http://www.adsciengineering.com/bpdcalc/index.php for example | |
362 | // http://www.adsciengineering.com/bpdcalc/index.php?n=20&p=0.25&start=0&stop=20&Submit=Generate | |
363 | // Appears to use at least 80-bit long double for 32 decimal digits accuracy, | |
364 | // but loses accuracy of display if leading zeros? | |
365 | // (if trailings zero then are exact values?) | |
366 | // so useful for testing 64-bit double accuracy. | |
367 | // P = 0.25, n = 20, k = 0 to 20 | |
368 | ||
369 | //0 C(20,0) * 0.25^0 * 0.75^20 0.00317121193893399322405457496643 | |
370 | //1 C(20,1) * 0.25^1 * 0.75^19 0.02114141292622662149369716644287 | |
371 | //2 C(20,2) * 0.25^2 * 0.75^18 0.06694780759971763473004102706909 | |
372 | //3 C(20,3) * 0.25^3 * 0.75^17 0.13389561519943526946008205413818 | |
373 | //4 C(20,4) * 0.25^4 * 0.75^16 0.18968545486586663173511624336242 | |
374 | //5 C(20,5) * 0.25^5 * 0.75^15 0.20233115185692440718412399291992 | |
375 | //6 C(20,6) * 0.25^6 * 0.75^14 0.16860929321410367265343666076660 | |
376 | //7 C(20,7) * 0.25^7 * 0.75^13 0.11240619547606911510229110717773 | |
377 | //8 C(20,8) * 0.25^8 * 0.75^12 0.06088668921620410401374101638793 | |
378 | //9 C(20,9) * 0.25^9 * 0.75^11 0.02706075076275737956166267395019 | |
379 | //10 C(20,10) * 0.25^10 * 0.75^10 0.00992227527967770583927631378173 | |
380 | //11 C(20,11) * 0.25^11 * 0.75^9 0.00300675008475081995129585266113 | |
381 | //12 C(20,12) * 0.25^12 * 0.75^8 0.00075168752118770498782396316528 | |
382 | //13 C(20,13) * 0.25^13 * 0.75^7 0.00015419231203850358724594116210 | |
383 | //14 C(20,14) * 0.25^14 * 0.75^6 0.00002569871867308393120765686035 | |
384 | //15 C(20,15) * 0.25^15 * 0.75^5 0.00000342649582307785749435424804 | |
385 | //16 C(20,16) * 0.25^16 * 0.75^4 0.00000035692664823727682232856750 | |
386 | //17 C(20,17) * 0.25^17 * 0.75^3 0.00000002799424692057073116302490 | |
387 | //18 C(20,18) * 0.25^18 * 0.75^2 0.00000000155523594003170728683471 | |
388 | //19 C(20,19) * 0.25^19 * 0.75^1 0.00000000005456968210637569427490 | |
389 | //20 C(20,20) * 0.25^20 * 0.75^0 0.00000000000090949470177292823791 | |
390 | ||
391 | ||
392 | BOOST_CHECK_CLOSE( | |
393 | pdf(binomial_distribution<RealType>(static_cast<RealType>(20), static_cast<RealType>(0.25)), | |
394 | static_cast<RealType>(10)), // k. | |
395 | static_cast<RealType>(0.00992227527967770583927631378173), // k=10 p = 0.25 | |
396 | tolerance); | |
397 | ||
398 | BOOST_CHECK_CLOSE( // k = 0 use different formula - only exp so more accurate. | |
399 | pdf(binomial_distribution<RealType>(static_cast<RealType>(20), static_cast<RealType>(0.25)), | |
400 | static_cast<RealType>(0)), // k. | |
401 | static_cast<RealType>(0.00317121193893399322405457496643), // k=0 p = 0.25 | |
402 | tolerance); | |
403 | ||
404 | BOOST_CHECK_CLOSE( // k = 20 use different formula - only exp so more accurate. | |
405 | pdf(binomial_distribution<RealType>(static_cast<RealType>(20), static_cast<RealType>(0.25)), | |
406 | static_cast<RealType>(20)), // k == n. | |
407 | static_cast<RealType>(0.00000000000090949470177292823791), // k=20 p = 0.25 | |
408 | tolerance); | |
409 | ||
410 | BOOST_CHECK_CLOSE( // k = 1. | |
411 | pdf(binomial_distribution<RealType>(static_cast<RealType>(20), static_cast<RealType>(0.25)), | |
412 | static_cast<RealType>(1)), // k. | |
413 | static_cast<RealType>(0.02114141292622662149369716644287), // k=1 p = 0.25 | |
414 | tolerance); | |
415 | ||
416 | // Some exact (probably) values. | |
417 | BOOST_CHECK_CLOSE( | |
418 | pdf(binomial_distribution<RealType>(static_cast<RealType>(8), static_cast<RealType>(0.25)), | |
419 | static_cast<RealType>(0)), // k. | |
420 | static_cast<RealType>(0.10011291503906250000000000000000), // k=0 p = 0.25 | |
421 | tolerance); | |
422 | ||
423 | BOOST_CHECK_CLOSE( // k = 1. | |
424 | pdf(binomial_distribution<RealType>(static_cast<RealType>(8), static_cast<RealType>(0.25)), | |
425 | static_cast<RealType>(1)), // k. | |
426 | static_cast<RealType>(0.26696777343750000000000000000000), // k=1 p = 0.25 | |
427 | tolerance); | |
428 | ||
429 | BOOST_CHECK_CLOSE( // k = 2. | |
430 | pdf(binomial_distribution<RealType>(static_cast<RealType>(8), static_cast<RealType>(0.25)), | |
431 | static_cast<RealType>(2)), // k. | |
432 | static_cast<RealType>(0.31146240234375000000000000000000), // k=2 p = 0.25 | |
433 | tolerance); | |
434 | ||
435 | BOOST_CHECK_CLOSE( // k = 3. | |
436 | pdf(binomial_distribution<RealType>(static_cast<RealType>(8), static_cast<RealType>(0.25)), | |
437 | static_cast<RealType>(3)), // k. | |
438 | static_cast<RealType>(0.20764160156250000000000000000000), // k=3 p = 0.25 | |
439 | tolerance); | |
440 | ||
441 | BOOST_CHECK_CLOSE( // k = 7. | |
442 | pdf(binomial_distribution<RealType>(static_cast<RealType>(8), static_cast<RealType>(0.25)), | |
443 | static_cast<RealType>(7)), // k. | |
444 | static_cast<RealType>(0.00036621093750000000000000000000), // k=7 p = 0.25 | |
445 | tolerance); | |
446 | ||
447 | BOOST_CHECK_CLOSE( // k = 8. | |
448 | pdf(binomial_distribution<RealType>(static_cast<RealType>(8), static_cast<RealType>(0.25)), | |
449 | static_cast<RealType>(8)), // k = n. | |
450 | static_cast<RealType>(0.00001525878906250000000000000000), // k=8 p = 0.25 | |
451 | tolerance); | |
452 | ||
453 | binomial_distribution<RealType> dist(static_cast<RealType>(8), static_cast<RealType>(0.25)); | |
454 | RealType x = static_cast<RealType>(0.125); | |
455 | using namespace std; // ADL of std names. | |
456 | // mean: | |
457 | BOOST_CHECK_CLOSE( | |
458 | mean(dist) | |
459 | , static_cast<RealType>(8 * 0.25), tol2); | |
460 | // variance: | |
461 | BOOST_CHECK_CLOSE( | |
462 | variance(dist) | |
463 | , static_cast<RealType>(8 * 0.25 * 0.75), tol2); | |
464 | // std deviation: | |
465 | BOOST_CHECK_CLOSE( | |
466 | standard_deviation(dist) | |
467 | , static_cast<RealType>(sqrt(8 * 0.25L * 0.75L)), tol2); | |
468 | // hazard: | |
469 | BOOST_CHECK_CLOSE( | |
470 | hazard(dist, x) | |
471 | , pdf(dist, x) / cdf(complement(dist, x)), tol2); | |
472 | // cumulative hazard: | |
473 | BOOST_CHECK_CLOSE( | |
474 | chf(dist, x) | |
475 | , -log(cdf(complement(dist, x))), tol2); | |
476 | // coefficient_of_variation: | |
477 | BOOST_CHECK_CLOSE( | |
478 | coefficient_of_variation(dist) | |
479 | , standard_deviation(dist) / mean(dist), tol2); | |
480 | // mode: | |
481 | BOOST_CHECK_CLOSE( | |
482 | mode(dist) | |
483 | , static_cast<RealType>(std::floor(9 * 0.25)), tol2); | |
484 | // skewness: | |
485 | BOOST_CHECK_CLOSE( | |
486 | skewness(dist) | |
487 | , static_cast<RealType>(0.40824829046386301636621401245098L), (std::max)(tol2, static_cast<RealType>(5e-29))); // test data has 32 digits only. | |
488 | // kurtosis: | |
489 | BOOST_CHECK_CLOSE( | |
490 | kurtosis(dist) | |
491 | , static_cast<RealType>(2.916666666666666666666666666666666666L), tol2); | |
492 | // kurtosis excess: | |
493 | BOOST_CHECK_CLOSE( | |
494 | kurtosis_excess(dist) | |
495 | , static_cast<RealType>(-0.08333333333333333333333333333333333333L), tol2); | |
496 | // Check kurtosis_excess == kurtosis -3; | |
497 | BOOST_CHECK_EQUAL(kurtosis(dist), static_cast<RealType>(3) + kurtosis_excess(dist)); | |
498 | ||
499 | // special cases for PDF: | |
500 | BOOST_CHECK_EQUAL( | |
501 | pdf( | |
502 | binomial_distribution<RealType>(static_cast<RealType>(8), static_cast<RealType>(0)), | |
503 | static_cast<RealType>(0)), static_cast<RealType>(1) | |
504 | ); | |
505 | BOOST_CHECK_EQUAL( | |
506 | pdf( | |
507 | binomial_distribution<RealType>(static_cast<RealType>(8), static_cast<RealType>(0)), | |
508 | static_cast<RealType>(0.0001)), static_cast<RealType>(0) | |
509 | ); | |
510 | BOOST_CHECK_EQUAL( | |
511 | pdf( | |
512 | binomial_distribution<RealType>(static_cast<RealType>(8), static_cast<RealType>(1)), | |
513 | static_cast<RealType>(0.001)), static_cast<RealType>(0) | |
514 | ); | |
515 | BOOST_CHECK_EQUAL( | |
516 | pdf( | |
517 | binomial_distribution<RealType>(static_cast<RealType>(8), static_cast<RealType>(1)), | |
518 | static_cast<RealType>(8)), static_cast<RealType>(1) | |
519 | ); | |
520 | BOOST_CHECK_EQUAL( | |
521 | pdf( | |
522 | binomial_distribution<RealType>(static_cast<RealType>(0), static_cast<RealType>(0.25)), | |
523 | static_cast<RealType>(0)), static_cast<RealType>(1) | |
524 | ); | |
525 | BOOST_MATH_CHECK_THROW( | |
526 | pdf( | |
527 | binomial_distribution<RealType>(static_cast<RealType>(-1), static_cast<RealType>(0.25)), | |
528 | static_cast<RealType>(0)), std::domain_error | |
529 | ); | |
530 | BOOST_MATH_CHECK_THROW( | |
531 | pdf( | |
532 | binomial_distribution<RealType>(static_cast<RealType>(8), static_cast<RealType>(-0.25)), | |
533 | static_cast<RealType>(0)), std::domain_error | |
534 | ); | |
535 | BOOST_MATH_CHECK_THROW( | |
536 | pdf( | |
537 | binomial_distribution<RealType>(static_cast<RealType>(8), static_cast<RealType>(1.25)), | |
538 | static_cast<RealType>(0)), std::domain_error | |
539 | ); | |
540 | BOOST_MATH_CHECK_THROW( | |
541 | pdf( | |
542 | binomial_distribution<RealType>(static_cast<RealType>(8), static_cast<RealType>(0.25)), | |
543 | static_cast<RealType>(-1)), std::domain_error | |
544 | ); | |
545 | BOOST_MATH_CHECK_THROW( | |
546 | pdf( | |
547 | binomial_distribution<RealType>(static_cast<RealType>(8), static_cast<RealType>(0.25)), | |
548 | static_cast<RealType>(9)), std::domain_error | |
549 | ); | |
550 | BOOST_MATH_CHECK_THROW( | |
551 | cdf( | |
552 | binomial_distribution<RealType>(static_cast<RealType>(8), static_cast<RealType>(0.25)), | |
553 | static_cast<RealType>(-1)), std::domain_error | |
554 | ); | |
555 | BOOST_MATH_CHECK_THROW( | |
556 | cdf( | |
557 | binomial_distribution<RealType>(static_cast<RealType>(8), static_cast<RealType>(0.25)), | |
558 | static_cast<RealType>(9)), std::domain_error | |
559 | ); | |
560 | BOOST_MATH_CHECK_THROW( | |
561 | cdf( | |
562 | binomial_distribution<RealType>(static_cast<RealType>(8), static_cast<RealType>(-0.25)), | |
563 | static_cast<RealType>(0)), std::domain_error | |
564 | ); | |
565 | BOOST_MATH_CHECK_THROW( | |
566 | cdf( | |
567 | binomial_distribution<RealType>(static_cast<RealType>(8), static_cast<RealType>(1.25)), | |
568 | static_cast<RealType>(0)), std::domain_error | |
569 | ); | |
570 | BOOST_MATH_CHECK_THROW( | |
571 | quantile( | |
572 | binomial_distribution<RealType>(static_cast<RealType>(8), static_cast<RealType>(-0.25)), | |
573 | static_cast<RealType>(0)), std::domain_error | |
574 | ); | |
575 | BOOST_MATH_CHECK_THROW( | |
576 | quantile( | |
577 | binomial_distribution<RealType>(static_cast<RealType>(8), static_cast<RealType>(1.25)), | |
578 | static_cast<RealType>(0)), std::domain_error | |
579 | ); | |
580 | ||
581 | BOOST_CHECK_EQUAL( | |
582 | quantile( | |
583 | binomial_distribution<RealType>(static_cast<RealType>(16), static_cast<RealType>(0.25)), | |
584 | static_cast<RealType>(0.01)), // Less than cdf == pdf(binomial_distribution<RealType>(16, 0.25), 0) | |
585 | static_cast<RealType>(0) // so expect zero as best approximation. | |
586 | ); | |
587 | ||
588 | BOOST_CHECK_EQUAL( | |
589 | cdf( | |
590 | binomial_distribution<RealType>(static_cast<RealType>(8), static_cast<RealType>(0.25)), | |
591 | static_cast<RealType>(8)), static_cast<RealType>(1) | |
592 | ); | |
593 | BOOST_CHECK_EQUAL( | |
594 | cdf( | |
595 | binomial_distribution<RealType>(static_cast<RealType>(8), static_cast<RealType>(0)), | |
596 | static_cast<RealType>(7)), static_cast<RealType>(1) | |
597 | ); | |
598 | BOOST_CHECK_EQUAL( | |
599 | cdf( | |
600 | binomial_distribution<RealType>(static_cast<RealType>(8), static_cast<RealType>(1)), | |
601 | static_cast<RealType>(7)), static_cast<RealType>(0) | |
602 | ); | |
603 | ||
604 | #endif | |
605 | ||
606 | { | |
607 | // This is a visual sanity check that everything is OK: | |
608 | binomial_distribution<RealType> my8dist(8., 0.25); // Note: double values (matching the distribution definition) avoid the need for any casting. | |
609 | //cout << "mean(my8dist) = " << boost::math::mean(my8dist) << endl; // mean(my8dist) = 2 | |
610 | //cout << "my8dist.trials() = " << my8dist.trials() << endl; // my8dist.trials() = 8 | |
611 | //cout << "my8dist.success_fraction() = " << my8dist.success_fraction() << endl; // my8dist.success_fraction() = 0.25 | |
612 | BOOST_CHECK_CLOSE(my8dist.trials(), static_cast<RealType>(8), tol2); | |
613 | BOOST_CHECK_CLOSE(my8dist.success_fraction(), static_cast<RealType>(0.25), tol2); | |
614 | ||
615 | //{ | |
616 | // int n = static_cast<int>(boost::math::tools::real_cast<double>(my8dist.trials())); | |
617 | // RealType sumcdf = 0.; | |
618 | // for (int k = 0; k <= n; k++) | |
619 | // { | |
620 | // cout << k << ' ' << pdf(my8dist, static_cast<RealType>(k)); | |
621 | // sumcdf += pdf(my8dist, static_cast<RealType>(k)); | |
622 | // cout << ' ' << sumcdf; | |
623 | // cout << ' ' << cdf(my8dist, static_cast<RealType>(k)); | |
624 | // cout << ' ' << sumcdf - cdf(my8dist, static_cast<RealType>(k)) << endl; | |
625 | // } // for k | |
626 | // } | |
627 | // n = 8, p =0.25 | |
628 | //k pdf cdf | |
629 | //0 0.1001129150390625 0.1001129150390625 | |
630 | //1 0.26696777343749994 0.36708068847656244 | |
631 | //2 0.31146240234375017 0.67854309082031261 | |
632 | //3 0.20764160156249989 0.8861846923828125 | |
633 | //4 0.086517333984375 0.9727020263671875 | |
634 | //5 0.023071289062499997 0.9957733154296875 | |
635 | //6 0.0038452148437500009 0.9996185302734375 | |
636 | //7 0.00036621093749999984 0.9999847412109375 | |
637 | //8 1.52587890625e-005 1 1 0 | |
638 | } | |
639 | #define T RealType | |
640 | #include "binomial_quantile.ipp" | |
641 | ||
642 | for(unsigned i = 0; i < binomial_quantile_data.size(); ++i) | |
643 | { | |
644 | using namespace boost::math::policies; | |
645 | typedef policy<discrete_quantile<boost::math::policies::real> > P1; | |
646 | typedef policy<discrete_quantile<integer_round_down> > P2; | |
647 | typedef policy<discrete_quantile<integer_round_up> > P3; | |
648 | typedef policy<discrete_quantile<integer_round_outwards> > P4; | |
649 | typedef policy<discrete_quantile<integer_round_inwards> > P5; | |
650 | typedef policy<discrete_quantile<integer_round_nearest> > P6; | |
651 | RealType tol = boost::math::tools::epsilon<RealType>() * 500; | |
652 | if(!boost::is_floating_point<RealType>::value) | |
653 | tol *= 10; // no lanczos approximation implies less accuracy | |
654 | RealType x; | |
655 | #if !defined(TEST_ROUNDING) || (TEST_ROUNDING == 1) | |
656 | // | |
657 | // Check full real value first: | |
658 | // | |
659 | binomial_distribution<RealType, P1> p1(binomial_quantile_data[i][0], binomial_quantile_data[i][1]); | |
660 | x = quantile(p1, binomial_quantile_data[i][2]); | |
661 | BOOST_CHECK_CLOSE_FRACTION(x, (RealType)binomial_quantile_data[i][3], tol); | |
662 | x = quantile(complement(p1, (RealType)binomial_quantile_data[i][2])); | |
663 | BOOST_CHECK_CLOSE_FRACTION(x, (RealType)binomial_quantile_data[i][4], tol); | |
664 | #endif | |
665 | #if !defined(TEST_ROUNDING) || (TEST_ROUNDING == 2) | |
666 | // | |
667 | // Now with round down to integer: | |
668 | // | |
669 | binomial_distribution<RealType, P2> p2(binomial_quantile_data[i][0], binomial_quantile_data[i][1]); | |
670 | x = quantile(p2, binomial_quantile_data[i][2]); | |
671 | BOOST_CHECK_EQUAL(x, (RealType)floor(binomial_quantile_data[i][3])); | |
672 | x = quantile(complement(p2, binomial_quantile_data[i][2])); | |
673 | BOOST_CHECK_EQUAL(x, (RealType)floor(binomial_quantile_data[i][4])); | |
674 | #endif | |
675 | #if !defined(TEST_ROUNDING) || (TEST_ROUNDING == 3) | |
676 | // | |
677 | // Now with round up to integer: | |
678 | // | |
679 | binomial_distribution<RealType, P3> p3(binomial_quantile_data[i][0], binomial_quantile_data[i][1]); | |
680 | x = quantile(p3, binomial_quantile_data[i][2]); | |
681 | BOOST_CHECK_EQUAL(x, (RealType)ceil(binomial_quantile_data[i][3])); | |
682 | x = quantile(complement(p3, binomial_quantile_data[i][2])); | |
683 | BOOST_CHECK_EQUAL(x, (RealType)ceil(binomial_quantile_data[i][4])); | |
684 | #endif | |
685 | #if !defined(TEST_ROUNDING) || (TEST_ROUNDING == 4) | |
686 | // | |
687 | // Now with round to integer "outside": | |
688 | // | |
689 | binomial_distribution<RealType, P4> p4(binomial_quantile_data[i][0], binomial_quantile_data[i][1]); | |
690 | x = quantile(p4, binomial_quantile_data[i][2]); | |
691 | BOOST_CHECK_EQUAL(x, (RealType)(binomial_quantile_data[i][2] < 0.5f ? floor(binomial_quantile_data[i][3]) : ceil(binomial_quantile_data[i][3]))); | |
692 | x = quantile(complement(p4, binomial_quantile_data[i][2])); | |
693 | BOOST_CHECK_EQUAL(x, (RealType)(binomial_quantile_data[i][2] < 0.5f ? ceil(binomial_quantile_data[i][4]) : floor(binomial_quantile_data[i][4]))); | |
694 | #endif | |
695 | #if !defined(TEST_ROUNDING) || (TEST_ROUNDING == 5) | |
696 | // | |
697 | // Now with round to integer "inside": | |
698 | // | |
699 | binomial_distribution<RealType, P5> p5(binomial_quantile_data[i][0], binomial_quantile_data[i][1]); | |
700 | x = quantile(p5, binomial_quantile_data[i][2]); | |
701 | BOOST_CHECK_EQUAL(x, (RealType)(binomial_quantile_data[i][2] < 0.5f ? ceil(binomial_quantile_data[i][3]) : floor(binomial_quantile_data[i][3]))); | |
702 | x = quantile(complement(p5, binomial_quantile_data[i][2])); | |
703 | BOOST_CHECK_EQUAL(x, (RealType)(binomial_quantile_data[i][2] < 0.5f ? floor(binomial_quantile_data[i][4]) : ceil(binomial_quantile_data[i][4]))); | |
704 | #endif | |
705 | #if !defined(TEST_ROUNDING) || (TEST_ROUNDING == 6) | |
706 | // | |
707 | // Now with round to nearest integer: | |
708 | // | |
709 | binomial_distribution<RealType, P6> p6(binomial_quantile_data[i][0], binomial_quantile_data[i][1]); | |
710 | x = quantile(p6, binomial_quantile_data[i][2]); | |
711 | BOOST_CHECK_EQUAL(x, (RealType)(floor(binomial_quantile_data[i][3] + 0.5f))); | |
712 | x = quantile(complement(p6, binomial_quantile_data[i][2])); | |
713 | BOOST_CHECK_EQUAL(x, (RealType)(floor(binomial_quantile_data[i][4] + 0.5f))); | |
714 | #endif | |
715 | } | |
716 | ||
717 | check_out_of_range<boost::math::binomial_distribution<RealType> >(1, 1); // (All) valid constructor parameter values. | |
718 | ||
719 | ||
720 | } // template <class RealType>void test_spots(RealType) | |
721 | ||
722 | BOOST_AUTO_TEST_CASE( test_main ) | |
723 | { | |
724 | BOOST_MATH_CONTROL_FP; | |
725 | // Check that can generate binomial distribution using one convenience methods: | |
726 | binomial_distribution<> mybn2(1., 0.5); // Using default RealType double. | |
727 | // but that | |
728 | // boost::math::binomial mybn1(1., 0.5); // Using typedef fails | |
729 | // error C2039: 'binomial' : is not a member of 'boost::math' | |
730 | ||
731 | // Basic sanity-check spot values. | |
732 | ||
733 | // (Parameter value, arbitrarily zero, only communicates the floating point type). | |
734 | #ifdef TEST_FLOAT | |
735 | test_spots(0.0F); // Test float. | |
736 | #endif | |
737 | #ifdef TEST_DOUBLE | |
738 | test_spots(0.0); // Test double. | |
739 | #endif | |
740 | #ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS | |
741 | #ifdef TEST_LDOUBLE | |
742 | test_spots(0.0L); // Test long double. | |
743 | #endif | |
744 | #if !defined(BOOST_MATH_NO_REAL_CONCEPT_TESTS) | |
745 | #ifdef TEST_REAL_CONCEPT | |
746 | test_spots(boost::math::concepts::real_concept(0.)); // Test real concept. | |
747 | #endif | |
748 | #endif | |
749 | #else | |
750 | std::cout << "<note>The long double tests have been disabled on this platform " | |
751 | "either because the long double overloads of the usual math functions are " | |
752 | "not available at all, or because they are too inaccurate for these tests " | |
753 | "to pass.</note>" << std::endl; | |
754 | #endif | |
755 | ||
756 | } // BOOST_AUTO_TEST_CASE( test_main ) | |
757 | ||
758 | /* | |
759 | ||
760 | Output is: | |
761 | ||
762 | Description: Autorun "J:\Cpp\MathToolkit\test\Math_test\Debug\test_binomial.exe" | |
763 | Running 1 test case... | |
764 | Tolerance for type float is 0.0119209 % | |
765 | Tolerance for type double is 2.22045e-011 % | |
766 | Tolerance for type long double is 2.22045e-011 % | |
767 | Tolerance for type class boost::math::concepts::real_concept is 2.22045e-011 % | |
768 | ||
769 | *** No errors detected | |
770 | ||
771 | ========== Build: 1 succeeded, 0 failed, 0 up-to-date, 0 skipped ========== | |
772 | ||
773 | ||
774 | */ |