Chris@16
|
1 // Copyright 2008 John Maddock
|
Chris@16
|
2 //
|
Chris@16
|
3 // Use, modification and distribution are subject to the
|
Chris@16
|
4 // Boost Software License, Version 1.0.
|
Chris@16
|
5 // (See accompanying file LICENSE_1_0.txt
|
Chris@16
|
6 // or copy at http://www.boost.org/LICENSE_1_0.txt)
|
Chris@16
|
7
|
Chris@16
|
8 #ifndef BOOST_MATH_DISTRIBUTIONS_DETAIL_HG_CDF_HPP
|
Chris@16
|
9 #define BOOST_MATH_DISTRIBUTIONS_DETAIL_HG_CDF_HPP
|
Chris@16
|
10
|
Chris@16
|
11 #include <boost/math/policies/error_handling.hpp>
|
Chris@16
|
12 #include <boost/math/distributions/detail/hypergeometric_pdf.hpp>
|
Chris@16
|
13
|
Chris@16
|
14 namespace boost{ namespace math{ namespace detail{
|
Chris@16
|
15
|
Chris@16
|
16 template <class T, class Policy>
|
Chris@16
|
17 T hypergeometric_cdf_imp(unsigned x, unsigned r, unsigned n, unsigned N, bool invert, const Policy& pol)
|
Chris@16
|
18 {
|
Chris@16
|
19 #ifdef BOOST_MSVC
|
Chris@16
|
20 # pragma warning(push)
|
Chris@16
|
21 # pragma warning(disable:4267)
|
Chris@16
|
22 #endif
|
Chris@16
|
23 BOOST_MATH_STD_USING
|
Chris@16
|
24 T result = 0;
|
Chris@16
|
25 T mode = floor(T(r + 1) * T(n + 1) / (N + 2));
|
Chris@16
|
26 if(x < mode)
|
Chris@16
|
27 {
|
Chris@16
|
28 result = hypergeometric_pdf<T>(x, r, n, N, pol);
|
Chris@16
|
29 T diff = result;
|
Chris@16
|
30 unsigned lower_limit = static_cast<unsigned>((std::max)(0, (int)(n + r) - (int)(N)));
|
Chris@16
|
31 while(diff > (invert ? T(1) : result) * tools::epsilon<T>())
|
Chris@16
|
32 {
|
Chris@16
|
33 diff = T(x) * T((N + x) - n - r) * diff / (T(1 + n - x) * T(1 + r - x));
|
Chris@16
|
34 result += diff;
|
Chris@16
|
35 BOOST_MATH_INSTRUMENT_VARIABLE(x);
|
Chris@16
|
36 BOOST_MATH_INSTRUMENT_VARIABLE(diff);
|
Chris@16
|
37 BOOST_MATH_INSTRUMENT_VARIABLE(result);
|
Chris@16
|
38 if(x == lower_limit)
|
Chris@16
|
39 break;
|
Chris@16
|
40 --x;
|
Chris@16
|
41 }
|
Chris@16
|
42 }
|
Chris@16
|
43 else
|
Chris@16
|
44 {
|
Chris@16
|
45 invert = !invert;
|
Chris@16
|
46 unsigned upper_limit = (std::min)(r, n);
|
Chris@16
|
47 if(x != upper_limit)
|
Chris@16
|
48 {
|
Chris@16
|
49 ++x;
|
Chris@16
|
50 result = hypergeometric_pdf<T>(x, r, n, N, pol);
|
Chris@16
|
51 T diff = result;
|
Chris@16
|
52 while((x <= upper_limit) && (diff > (invert ? T(1) : result) * tools::epsilon<T>()))
|
Chris@16
|
53 {
|
Chris@16
|
54 diff = T(n - x) * T(r - x) * diff / (T(x + 1) * T((N + x + 1) - n - r));
|
Chris@16
|
55 result += diff;
|
Chris@16
|
56 ++x;
|
Chris@16
|
57 BOOST_MATH_INSTRUMENT_VARIABLE(x);
|
Chris@16
|
58 BOOST_MATH_INSTRUMENT_VARIABLE(diff);
|
Chris@16
|
59 BOOST_MATH_INSTRUMENT_VARIABLE(result);
|
Chris@16
|
60 }
|
Chris@16
|
61 }
|
Chris@16
|
62 }
|
Chris@16
|
63 if(invert)
|
Chris@16
|
64 result = 1 - result;
|
Chris@16
|
65 return result;
|
Chris@16
|
66 #ifdef BOOST_MSVC
|
Chris@16
|
67 # pragma warning(pop)
|
Chris@16
|
68 #endif
|
Chris@16
|
69 }
|
Chris@16
|
70
|
Chris@16
|
71 template <class T, class Policy>
|
Chris@16
|
72 inline T hypergeometric_cdf(unsigned x, unsigned r, unsigned n, unsigned N, bool invert, const Policy&)
|
Chris@16
|
73 {
|
Chris@16
|
74 BOOST_FPU_EXCEPTION_GUARD
|
Chris@16
|
75 typedef typename tools::promote_args<T>::type result_type;
|
Chris@16
|
76 typedef typename policies::evaluation<result_type, Policy>::type value_type;
|
Chris@16
|
77 typedef typename policies::normalise<
|
Chris@16
|
78 Policy,
|
Chris@16
|
79 policies::promote_float<false>,
|
Chris@16
|
80 policies::promote_double<false>,
|
Chris@16
|
81 policies::discrete_quantile<>,
|
Chris@16
|
82 policies::assert_undefined<> >::type forwarding_policy;
|
Chris@16
|
83
|
Chris@16
|
84 value_type result;
|
Chris@16
|
85 result = detail::hypergeometric_cdf_imp<value_type>(x, r, n, N, invert, forwarding_policy());
|
Chris@16
|
86 if(result > 1)
|
Chris@16
|
87 {
|
Chris@16
|
88 result = 1;
|
Chris@16
|
89 }
|
Chris@16
|
90 if(result < 0)
|
Chris@16
|
91 {
|
Chris@16
|
92 result = 0;
|
Chris@16
|
93 }
|
Chris@16
|
94 return policies::checked_narrowing_cast<result_type, forwarding_policy>(result, "boost::math::hypergeometric_cdf<%1%>(%1%,%1%,%1%,%1%)");
|
Chris@16
|
95 }
|
Chris@16
|
96
|
Chris@16
|
97 }}} // namespaces
|
Chris@16
|
98
|
Chris@16
|
99 #endif
|
Chris@16
|
100
|