1 // Copyright John Maddock 2008
2 // Copyright Paul A. Bristow
3 // Copyright Gautam Sewani
5 // Use, modification and distribution are subject to the
6 // Boost Software License, Version 1.0.
7 // (See accompanying file LICENSE_1_0.txt
8 // or copy at http://www.boost.org/LICENSE_1_0.txt)
10 #define BOOST_MATH_OVERFLOW_ERROR_POLICY throw_on_error
11 #include <boost/math/concepts/real_concept.hpp> // for real_concept
12 #include <boost/math/distributions/hypergeometric.hpp>
14 #define BOOST_TEST_MAIN
15 #include <boost/test/unit_test.hpp> // Boost.Test
16 #include <boost/test/results_collector.hpp>
17 #include <boost/test/unit_test.hpp>
18 #include <boost/test/tools/floating_point_comparison.hpp>
23 using std::setprecision
;
25 #include <boost/array.hpp>
26 #include "functor.hpp"
27 #include "handle_test_result.hpp"
28 #include "table_type.hpp"
30 #define BOOST_CHECK_EX(a) \
32 unsigned int failures = boost::unit_test::results_collector.results( boost::unit_test::framework::current_test_case().p_id ).p_assertions_failed;\
34 if(failures != boost::unit_test::results_collector.results( boost::unit_test::framework::current_test_case().p_id ).p_assertions_failed)\
36 std::cerr << "Failure was with data ";\
37 std::cerr << std::setprecision(35); \
38 std::cerr << "x = " << x << ", r = " << r << ", n = " << n\
39 << ", N = " << N << ", p = " << cp << ", q = " << ccp << std::endl;\
43 void expected_results()
46 // Define the max and mean errors expected for
47 // various compilers and platforms.
49 const char* largest_type
;
50 #ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
51 if(boost::math::policies::digits
<double, boost::math::policies::policy
<> >() == boost::math::policies::digits
<long double, boost::math::policies::policy
<> >())
53 largest_type
= "(long\\s+)?double|real_concept";
57 largest_type
= "long double|real_concept";
60 largest_type
= "(long\\s+)?double";
62 #ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
63 if((boost::math::tools::digits
<long double>() > boost::math::tools::digits
<double>())
64 && (boost::math::tools::digits
<long double>() < 100))
67 // Some split of errors from long double into double:
73 "double", // test type(s)
74 "Random.*", // test data group
75 ".*", 1500, 1500); // test function
80 "double", // test type(s)
81 ".*", // test data group
82 ".*", 10, 10); // test function
90 "real_concept", // test type(s)
91 "Random.*", // test data group
92 ".*", 250000000, 25000000); // test function
97 largest_type
, // test type(s)
98 "Random.*", // test data group
99 ".*", 10000000, 5000000); // test function
104 largest_type
, // test type(s)
105 ".*", // test data group
106 ".*", 50, 20); // test function
110 inline unsigned make_unsigned(T x
)
112 return static_cast<unsigned>(x
);
115 inline unsigned make_unsigned(boost::math::concepts::real_concept x
)
117 return static_cast<unsigned>(x
.value());
121 T
pdf_tester(T r
, T n
, T N
, T x
)
123 boost::math::hypergeometric_distribution
<T
> d(make_unsigned(r
), make_unsigned(n
), make_unsigned(N
));
128 T
cdf_tester(T r
, T n
, T N
, T x
)
130 boost::math::hypergeometric_distribution
<T
> d(make_unsigned(r
), make_unsigned(n
), make_unsigned(N
));
135 T
ccdf_tester(T r
, T n
, T N
, T x
)
137 boost::math::hypergeometric_distribution
<T
> d(make_unsigned(r
), make_unsigned(n
), make_unsigned(N
));
138 return cdf(complement(d
, x
));
141 template <class Real
, class T
>
142 void do_test_hypergeometric(const T
& data
, const char* type_name
, const char* test_name
)
144 // warning suppression:
149 #if !defined(TEST_QUANT) || (TEST_QUANT == 0)
150 typedef Real value_type
;
152 typedef value_type (*pg
)(value_type
, value_type
, value_type
, value_type
);
153 #if defined(BOOST_MATH_NO_DEDUCED_FUNCTION_POINTERS)
154 pg funcp
= pdf_tester
<value_type
>;
156 pg funcp
= pdf_tester
;
159 boost::math::tools::test_result
<value_type
> result
;
161 std::cout
<< "Testing " << test_name
<< " with type " << type_name
162 << "\n~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\n";
165 // test hypergeometric against data:
167 result
= boost::math::tools::test_hetero
<Real
>(
169 bind_func
<Real
>(funcp
, 0, 1, 2, 3),
170 extract_result
<Real
>(4));
171 handle_test_result(result
, data
[result
.worst()], result
.worst(), type_name
, "hypergeometric PDF", test_name
);
173 #if defined(BOOST_MATH_NO_DEDUCED_FUNCTION_POINTERS)
174 funcp
= cdf_tester
<value_type
>;
180 // test hypergeometric against data:
182 result
= boost::math::tools::test_hetero
<Real
>(
184 bind_func
<Real
>(funcp
, 0, 1, 2, 3),
185 extract_result
<Real
>(5));
186 handle_test_result(result
, data
[result
.worst()], result
.worst(), type_name
, "hypergeometric CDF", test_name
);
188 #if defined(BOOST_MATH_NO_DEDUCED_FUNCTION_POINTERS)
189 funcp
= ccdf_tester
<value_type
>;
195 // test hypergeometric against data:
197 result
= boost::math::tools::test_hetero
<Real
>(
199 bind_func
<Real
>(funcp
, 0, 1, 2, 3),
200 extract_result
<Real
>(6));
201 handle_test_result(result
, data
[result
.worst()], result
.worst(), type_name
, "hypergeometric CDF complement", test_name
);
202 std::cout
<< std::endl
;
206 template <class Real
, class T
>
207 void do_test_hypergeometric_quantile(const T
& data
, const char* type_name
, const char* test_name
)
209 typedef Real value_type
;
211 std::cout
<< "Checking quantiles with " << test_name
<< " with type " << type_name
212 << "\n~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\n";
214 if(boost::math::tools::digits
<value_type
>() > 50)
216 for(unsigned i
= 0; i
< data
.size(); ++i
)
218 using namespace boost::math::policies
;
220 unsigned r
= make_unsigned(data
[i
][0]);
221 unsigned n
= make_unsigned(data
[i
][1]);
222 unsigned N
= make_unsigned(data
[i
][2]);
223 unsigned x
= make_unsigned(data
[i
][3]);
224 value_type cp
= data
[i
][5];
225 value_type ccp
= data
[i
][6];
227 // A bit of warning suppression:
236 #if !defined(TEST_QUANT) || (TEST_QUANT == 1)
237 boost::math::hypergeometric_distribution
<value_type
,
238 policy
<discrete_quantile
<integer_round_up
> > > du(r
, n
, N
);
240 if((cp
< 0.9) && (cp
> boost::math::tools::min_value
<value_type
>()))
242 BOOST_CHECK_EX(quantile(du
, cp
) >= x
);
244 if((ccp
< 0.9) && (ccp
> boost::math::tools::min_value
<value_type
>()))
246 BOOST_CHECK_EX(quantile(complement(du
, ccp
)) >= x
);
249 #if !defined(TEST_QUANT) || (TEST_QUANT == 2)
250 boost::math::hypergeometric_distribution
<value_type
,
251 policy
<discrete_quantile
<integer_round_down
> > > dl(r
, n
, N
);
252 if((cp
< 0.9) && (cp
> boost::math::tools::min_value
<value_type
>()))
254 BOOST_CHECK_EX(quantile(dl
, cp
) <= x
);
256 if((ccp
< 0.9) && (ccp
> boost::math::tools::min_value
<value_type
>()))
258 BOOST_CHECK_EX(quantile(complement(dl
, ccp
)) <= x
);
261 #if !defined(TEST_QUANT) || (TEST_QUANT == 3)
262 boost::math::hypergeometric_distribution
<value_type
,
263 policy
<discrete_quantile
<integer_round_nearest
> > > dn(r
, n
, N
);
265 if((cp
< 0.9) && (cp
> boost::math::tools::min_value
<value_type
>()))
267 BOOST_CHECK_EX(quantile(dn
, cp
) == x
);
269 if((ccp
< 0.9) && (ccp
> boost::math::tools::min_value
<value_type
>()))
271 BOOST_CHECK_EX(quantile(complement(dn
, ccp
)) == x
);
274 #if !defined(TEST_QUANT) || (TEST_QUANT == 4)
275 boost::math::hypergeometric_distribution
<value_type
,
276 policy
<discrete_quantile
<integer_round_outwards
> > > dou(r
, n
, N
);
278 if((cp
< 0.9) && (cp
> boost::math::tools::min_value
<value_type
>()))
282 BOOST_CHECK_EX(quantile(dou
, cp
) <= x
);
286 BOOST_CHECK_EX(quantile(dou
, cp
) >= x
);
289 if((ccp
< 0.9) && (ccp
> boost::math::tools::min_value
<value_type
>()))
293 BOOST_CHECK_EX(quantile(complement(dou
, ccp
)) >= x
);
297 BOOST_CHECK_EX(quantile(complement(dou
, ccp
)) <= x
);
301 #if !defined(TEST_QUANT) || (TEST_QUANT == 5)
302 boost::math::hypergeometric_distribution
<value_type
,
303 policy
<discrete_quantile
<integer_round_inwards
> > > di(r
, n
, N
);
305 if((cp
< 0.9) && (cp
> boost::math::tools::min_value
<value_type
>()))
309 BOOST_CHECK_EX(quantile(di
, cp
) >= x
);
313 BOOST_CHECK_EX(quantile(di
, cp
) <= x
);
316 if((ccp
< 0.9) && (ccp
> boost::math::tools::min_value
<value_type
>()))
320 BOOST_CHECK_EX(quantile(complement(di
, ccp
)) <= x
);
324 BOOST_CHECK_EX(quantile(complement(di
, ccp
)) >= x
);
333 template <class RealType
>
334 void test_spot(unsigned x
, unsigned n
, unsigned r
, unsigned N
,
335 RealType p
, RealType cp
, RealType ccp
, RealType tol
)
338 // A bit of warning suppression:
348 #if !defined(TEST_QUANT) || (TEST_QUANT == 0)
349 boost::math::hypergeometric_distribution
<RealType
> d(r
, n
, N
);
351 std::pair
<unsigned, unsigned> extent
= range(d
);
353 BOOST_CHECK_CLOSE(pdf(d
, x
), p
, tol
);
354 BOOST_CHECK_CLOSE(cdf(d
, x
), cp
, tol
);
355 BOOST_CHECK_CLOSE(cdf(complement(d
, x
)), ccp
, tol
);
356 // Again with real-value arguments:
357 BOOST_CHECK_CLOSE(pdf(d
, static_cast<RealType
>(x
)), p
, tol
);
358 BOOST_CHECK_CLOSE(cdf(d
, static_cast<RealType
>(x
)), cp
, tol
);
359 BOOST_CHECK_CLOSE(cdf(complement(d
, static_cast<RealType
>(x
))), ccp
, tol
);
361 // Quantiles, don't bother checking these for type float
362 // as there's not enough precision in the p and q values
363 // to get back to where we started:
365 if(boost::math::tools::digits
<RealType
>() > 50)
367 using namespace boost::math::policies
;
369 boost::math::hypergeometric_distribution
<RealType
,
370 policy
<discrete_quantile
<integer_round_up
> > > du(r
, n
, N
);
372 BOOST_CHECK_EX(quantile(du
, cp
) >= x
);
373 BOOST_CHECK_EX(quantile(complement(du
, ccp
)) >= x
);
375 boost::math::hypergeometric_distribution
<RealType
,
376 policy
<discrete_quantile
<integer_round_down
> > > dl(r
, n
, N
);
377 BOOST_CHECK_EX(quantile(dl
, cp
) <= x
);
378 BOOST_CHECK_EX(quantile(complement(dl
, ccp
)) <= x
);
380 boost::math::hypergeometric_distribution
<RealType
,
381 policy
<discrete_quantile
<integer_round_nearest
> > > dn(r
, n
, N
);
383 BOOST_CHECK_EX(quantile(dn
, cp
) == x
);
384 BOOST_CHECK_EX(quantile(complement(dn
, ccp
)) == x
);
387 // Error checking of out of bounds arguments:
389 BOOST_MATH_CHECK_THROW(pdf(d
, extent
.second
+ 1), std::domain_error
);
390 BOOST_MATH_CHECK_THROW(cdf(d
, extent
.second
+ 1), std::domain_error
);
391 BOOST_MATH_CHECK_THROW(cdf(complement(d
, extent
.second
+ 1)), std::domain_error
);
394 BOOST_MATH_CHECK_THROW(pdf(d
, extent
.first
- 1), std::domain_error
);
395 BOOST_MATH_CHECK_THROW(cdf(d
, extent
.first
- 1), std::domain_error
);
396 BOOST_MATH_CHECK_THROW(cdf(complement(d
, extent
.first
- 1)), std::domain_error
);
398 BOOST_MATH_CHECK_THROW(quantile(d
, 1.1f
), std::domain_error
);
399 BOOST_MATH_CHECK_THROW(quantile(complement(d
, 1.1f
)), std::domain_error
);
400 BOOST_MATH_CHECK_THROW(quantile(d
, -0.001f
), std::domain_error
);
401 BOOST_MATH_CHECK_THROW(quantile(complement(d
, -0.001f
)), std::domain_error
);
403 // Checking of extreme values:
405 BOOST_CHECK_EQUAL(quantile(d
, 0), extent
.first
);
406 BOOST_CHECK_EQUAL(quantile(d
, 1), extent
.second
);
407 BOOST_CHECK_EQUAL(quantile(complement(d
, 0)), extent
.second
);
408 BOOST_CHECK_EQUAL(quantile(complement(d
, 1)), extent
.first
);
409 BOOST_CHECK_EQUAL(cdf(d
, extent
.second
), 1);
410 BOOST_CHECK_EQUAL(cdf(complement(d
, extent
.second
)), 0);
414 template <class RealType
>
415 void test_spots(RealType
/*T*/, const char* type_name
)
417 // Basic sanity checks.
418 // 50 eps as a percentage, up to a maximum of double precision
419 // Test data taken from Mathematica 6
421 #include "hypergeometric_test_data.ipp"
422 do_test_hypergeometric
<T
>(hypergeometric_test_data
, type_name
, "Mathematica data");
424 #include "hypergeometric_dist_data2.ipp"
425 if(boost::is_floating_point
<RealType
>::value
)
428 // Don't test this for real_concept: it's too slow!!!
430 do_test_hypergeometric
<T
>(hypergeometric_dist_data2
, type_name
, "Random large data");
433 do_test_hypergeometric_quantile
<T
>(hypergeometric_test_data
, type_name
, "Mathematica data");
434 if(boost::is_floating_point
<RealType
>::value
)
437 // Don't test this for real_concept: it's too slow!!!
439 do_test_hypergeometric_quantile
<T
>(hypergeometric_dist_data2
, type_name
, "Random large data");
442 RealType tolerance
= (std::max
)(
443 static_cast<RealType
>(2e-16L), // limit of test data
444 boost::math::tools::epsilon
<RealType
>());
445 cout
<<"Absolute tolerance:"<<tolerance
<<endl
;
447 tolerance
*= 50 * 100; // 50eps as a persentage
448 cout
<< "Tolerance for type " << typeid(RealType
).name() << " is " << tolerance
<< " %" << endl
;
451 // These sanity check values were calculated using the online
452 // calculator at http://stattrek.com/Tables/Hypergeometric.aspx
453 // It's assumed that the test values are accurate to no more than
456 test_spot(20, 200, 50, 500, static_cast<T
>(0.120748236361163), static_cast<T
>(0.563532430195156), static_cast<T
>(1 - 0.563532430195156), tolerance
);
457 test_spot(53, 452, 64, 500, static_cast<T
>(0.0184749573044286), static_cast<T
>(0.0299118078796907), static_cast<T
>(1 - 0.0299118078796907), tolerance
);
458 test_spot(32, 1287, 128, 5000, static_cast<T
>(0.0807012167418264), static_cast<T
>(0.469768774237742), static_cast<T
>(1 - 0.469768774237742), tolerance
);
459 test_spot(1, 13, 4, 26, static_cast<T
>(0.248695652173913), static_cast<T
>(0.296521739130435), static_cast<T
>(1 - 0.296521739130435), tolerance
);
460 test_spot(2, 13, 4, 26, static_cast<T
>(0.40695652173913), static_cast<T
>(0.703478260869565), static_cast<T
>(1 - 0.703478260869565), tolerance
);
461 test_spot(3, 13, 4, 26, static_cast<T
>(0.248695652173913), static_cast<T
>(0.952173913043478), static_cast<T
>(1 - 0.952173913043478), tolerance
);
462 test_spot(40, 70, 89, 170, static_cast<T
>(0.0721901023798991), static_cast<T
>(0.885447799131944), static_cast<T
>(1 - 0.885447799131944), tolerance
);
464 boost::math::hypergeometric_distribution
<RealType
> d(50, 200, 500);
465 BOOST_CHECK_EQUAL(range(d
).first
, 0u);
466 BOOST_CHECK_EQUAL(range(d
).second
, 50u);
467 BOOST_CHECK_CLOSE(mean(d
), static_cast<RealType
>(20), tolerance
);
468 BOOST_CHECK_CLOSE(mode(d
), static_cast<RealType
>(20), tolerance
);
469 BOOST_CHECK_CLOSE(variance(d
), static_cast<RealType
>(10.821643286573146292585170340681L), tolerance
);
470 BOOST_CHECK_CLOSE(skewness(d
), static_cast<RealType
>(0.048833071022952084732902910189366L), tolerance
);
471 BOOST_CHECK_CLOSE(kurtosis_excess(d
), static_cast<RealType
>(2.5155486690782804816404001878293L), tolerance
);
472 BOOST_CHECK_CLOSE(kurtosis(d
), kurtosis_excess(d
) + 3, tolerance
);
473 BOOST_CHECK_EQUAL(quantile(d
, 0.5f
), median(d
));
475 BOOST_MATH_CHECK_THROW(d
= boost::math::hypergeometric_distribution
<RealType
>(501, 40, 500), std::domain_error
);
476 BOOST_MATH_CHECK_THROW(d
= boost::math::hypergeometric_distribution
<RealType
>(40, 501, 500), std::domain_error
);
480 BOOST_AUTO_TEST_CASE( test_main
)
483 // Basic sanity-check spot values.
484 // (Parameter value, arbitrarily zero, only communicates the floating point type).
485 test_spots(0.0F
, "float"); // Test float. OK at decdigits = 0 tolerance = 0.0001 %
486 test_spots(0.0, "double"); // Test double. OK at decdigits 7, tolerance = 1e07 %
487 #ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
488 test_spots(0.0L, "long double"); // Test long double.
489 #ifndef BOOST_MATH_NO_REAL_CONCEPT_TESTS
490 test_spots(boost::math::concepts::real_concept(0), "real_concept"); // Test real_concept.
493 std::cout
<< "<note>The long double tests have been disabled on this platform "
494 "either because the long double overloads of the usual math functions are "
495 "not available at all, or because they are too inaccurate for these tests "
496 "to pass.</note>" << std::endl
;
500 } // BOOST_AUTO_TEST_CASE( test_main )