| // Copyright 2008 John Maddock |
| // |
| // Use, modification and distribution are subject to the |
| // Boost Software License, Version 1.0. |
| // (See accompanying file LICENSE_1_0.txt |
| // or copy at http://www.boost.org/LICENSE_1_0.txt) |
| |
| #ifndef BOOST_MATH_DISTRIBUTIONS_DETAIL_HG_QUANTILE_HPP |
| #define BOOST_MATH_DISTRIBUTIONS_DETAIL_HG_QUANTILE_HPP |
| |
| #include <boost/math/policies/error_handling.hpp> |
| #include <boost/math/distributions/detail/hypergeometric_pdf.hpp> |
| |
| namespace boost{ namespace math{ namespace detail{ |
| |
| template <class T> |
| inline unsigned round_x_from_p(unsigned x, T p, T cum, T fudge_factor, unsigned lbound, unsigned /*ubound*/, const policies::discrete_quantile<policies::integer_round_down>&) |
| { |
| if((p < cum * fudge_factor) && (x != lbound)) |
| { |
| BOOST_MATH_INSTRUMENT_VARIABLE(x-1); |
| return --x; |
| } |
| return x; |
| } |
| |
| template <class T> |
| inline unsigned round_x_from_p(unsigned x, T p, T cum, T fudge_factor, unsigned /*lbound*/, unsigned ubound, const policies::discrete_quantile<policies::integer_round_up>&) |
| { |
| if((cum < p * fudge_factor) && (x != ubound)) |
| { |
| BOOST_MATH_INSTRUMENT_VARIABLE(x+1); |
| return ++x; |
| } |
| return x; |
| } |
| |
| template <class T> |
| inline unsigned round_x_from_p(unsigned x, T p, T cum, T fudge_factor, unsigned lbound, unsigned ubound, const policies::discrete_quantile<policies::integer_round_inwards>&) |
| { |
| if(p >= 0.5) |
| return round_x_from_p(x, p, cum, fudge_factor, lbound, ubound, policies::discrete_quantile<policies::integer_round_down>()); |
| return round_x_from_p(x, p, cum, fudge_factor, lbound, ubound, policies::discrete_quantile<policies::integer_round_up>()); |
| } |
| |
| template <class T> |
| inline unsigned round_x_from_p(unsigned x, T p, T cum, T fudge_factor, unsigned lbound, unsigned ubound, const policies::discrete_quantile<policies::integer_round_outwards>&) |
| { |
| if(p >= 0.5) |
| return round_x_from_p(x, p, cum, fudge_factor, lbound, ubound, policies::discrete_quantile<policies::integer_round_up>()); |
| return round_x_from_p(x, p, cum, fudge_factor, lbound, ubound, policies::discrete_quantile<policies::integer_round_down>()); |
| } |
| |
| template <class T> |
| inline unsigned round_x_from_p(unsigned x, T /*p*/, T /*cum*/, T /*fudge_factor*/, unsigned /*lbound*/, unsigned /*ubound*/, const policies::discrete_quantile<policies::integer_round_nearest>&) |
| { |
| return x; |
| } |
| |
| template <class T> |
| inline unsigned round_x_from_q(unsigned x, T q, T cum, T fudge_factor, unsigned lbound, unsigned /*ubound*/, const policies::discrete_quantile<policies::integer_round_down>&) |
| { |
| if((q * fudge_factor > cum) && (x != lbound)) |
| { |
| BOOST_MATH_INSTRUMENT_VARIABLE(x-1); |
| return --x; |
| } |
| return x; |
| } |
| |
| template <class T> |
| inline unsigned round_x_from_q(unsigned x, T q, T cum, T fudge_factor, unsigned /*lbound*/, unsigned ubound, const policies::discrete_quantile<policies::integer_round_up>&) |
| { |
| if((q < cum * fudge_factor) && (x != ubound)) |
| { |
| BOOST_MATH_INSTRUMENT_VARIABLE(x+1); |
| return ++x; |
| } |
| return x; |
| } |
| |
| template <class T> |
| inline unsigned round_x_from_q(unsigned x, T q, T cum, T fudge_factor, unsigned lbound, unsigned ubound, const policies::discrete_quantile<policies::integer_round_inwards>&) |
| { |
| if(q < 0.5) |
| return round_x_from_q(x, q, cum, fudge_factor, lbound, ubound, policies::discrete_quantile<policies::integer_round_down>()); |
| return round_x_from_q(x, q, cum, fudge_factor, lbound, ubound, policies::discrete_quantile<policies::integer_round_up>()); |
| } |
| |
| template <class T> |
| inline unsigned round_x_from_q(unsigned x, T q, T cum, T fudge_factor, unsigned lbound, unsigned ubound, const policies::discrete_quantile<policies::integer_round_outwards>&) |
| { |
| if(q >= 0.5) |
| return round_x_from_q(x, q, cum, fudge_factor, lbound, ubound, policies::discrete_quantile<policies::integer_round_down>()); |
| return round_x_from_q(x, q, cum, fudge_factor, lbound, ubound, policies::discrete_quantile<policies::integer_round_up>()); |
| } |
| |
| template <class T> |
| inline unsigned round_x_from_q(unsigned x, T /*q*/, T /*cum*/, T /*fudge_factor*/, unsigned /*lbound*/, unsigned /*ubound*/, const policies::discrete_quantile<policies::integer_round_nearest>&) |
| { |
| return x; |
| } |
| |
| template <class T, class Policy> |
| unsigned hypergeometric_quantile_imp(T p, T q, unsigned r, unsigned n, unsigned N, const Policy& pol) |
| { |
| #ifdef BOOST_MSVC |
| # pragma warning(push) |
| # pragma warning(disable:4267) |
| #endif |
| typedef typename Policy::discrete_quantile_type discrete_quantile_type; |
| BOOST_MATH_STD_USING |
| BOOST_FPU_EXCEPTION_GUARD |
| T result; |
| T fudge_factor = 1 + tools::epsilon<T>() * ((N <= boost::math::prime(boost::math::max_prime - 1)) ? 50 : 2 * N); |
| unsigned base = static_cast<unsigned>((std::max)(0, (int)(n + r) - (int)(N))); |
| unsigned lim = (std::min)(r, n); |
| |
| BOOST_MATH_INSTRUMENT_VARIABLE(p); |
| BOOST_MATH_INSTRUMENT_VARIABLE(q); |
| BOOST_MATH_INSTRUMENT_VARIABLE(r); |
| BOOST_MATH_INSTRUMENT_VARIABLE(n); |
| BOOST_MATH_INSTRUMENT_VARIABLE(N); |
| BOOST_MATH_INSTRUMENT_VARIABLE(fudge_factor); |
| BOOST_MATH_INSTRUMENT_VARIABLE(base); |
| BOOST_MATH_INSTRUMENT_VARIABLE(lim); |
| |
| if(p <= 0.5) |
| { |
| unsigned x = base; |
| result = hypergeometric_pdf<T>(x, r, n, N, pol); |
| T diff = result; |
| while(result < p) |
| { |
| diff = (diff > tools::min_value<T>() * 8) |
| ? T(n - x) * T(r - x) * diff / (T(x + 1) * T(N + x + 1 - n - r)) |
| : hypergeometric_pdf<T>(x + 1, r, n, N, pol); |
| if(result + diff / 2 > p) |
| break; |
| ++x; |
| result += diff; |
| #ifdef BOOST_MATH_INSTRUMENT |
| if(diff != 0) |
| { |
| BOOST_MATH_INSTRUMENT_VARIABLE(x); |
| BOOST_MATH_INSTRUMENT_VARIABLE(diff); |
| BOOST_MATH_INSTRUMENT_VARIABLE(result); |
| } |
| #endif |
| } |
| return round_x_from_p(x, p, result, fudge_factor, base, lim, discrete_quantile_type()); |
| } |
| else |
| { |
| unsigned x = lim; |
| result = 0; |
| T diff = hypergeometric_pdf<T>(x, r, n, N, pol); |
| while(result + diff / 2 < q) |
| { |
| result += diff; |
| diff = (diff > tools::min_value<T>() * 8) |
| ? x * T(N + x - n - r) * diff / (T(1 + n - x) * T(1 + r - x)) |
| : hypergeometric_pdf<T>(x - 1, r, n, N, pol); |
| --x; |
| #ifdef BOOST_MATH_INSTRUMENT |
| if(diff != 0) |
| { |
| BOOST_MATH_INSTRUMENT_VARIABLE(x); |
| BOOST_MATH_INSTRUMENT_VARIABLE(diff); |
| BOOST_MATH_INSTRUMENT_VARIABLE(result); |
| } |
| #endif |
| } |
| return round_x_from_q(x, q, result, fudge_factor, base, lim, discrete_quantile_type()); |
| } |
| #ifdef BOOST_MSVC |
| # pragma warning(pop) |
| #endif |
| } |
| |
| template <class T, class Policy> |
| inline unsigned hypergeometric_quantile(T p, T q, unsigned r, unsigned n, unsigned N, const Policy&) |
| { |
| BOOST_FPU_EXCEPTION_GUARD |
| typedef typename tools::promote_args<T>::type result_type; |
| typedef typename policies::evaluation<result_type, Policy>::type value_type; |
| typedef typename policies::normalise< |
| Policy, |
| policies::promote_float<false>, |
| policies::promote_double<false>, |
| policies::assert_undefined<> >::type forwarding_policy; |
| |
| return detail::hypergeometric_quantile_imp<value_type>(p, q, r, n, N, forwarding_policy()); |
| } |
| |
| }}} // namespaces |
| |
| #endif |
| |