distributions/detail/hypergeometric_quantile.hpp

0001 // Copyright 2008 John Maddock
0002 //
0003 // Use, modification and distribution are subject to the
0004 // Boost Software License, Version 1.0.
0005 // (See accompanying file LICENSE_1_0.txt
0006 // or copy at http://www.boost.org/LICENSE_1_0.txt)
0007
0008 #ifndef BOOST_MATH_DISTRIBUTIONS_DETAIL_HG_QUANTILE_HPP
0009 #define BOOST_MATH_DISTRIBUTIONS_DETAIL_HG_QUANTILE_HPP
0010
0011 #include <boost/math/policies/error_handling.hpp>
0012 #include <boost/math/distributions/detail/hypergeometric_pdf.hpp>
0013
0014 namespace boost{ namespace math{ namespace detail{
0015
0016 template <class T>
0017 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>&)
0018 {
0019    if((p < cum * fudge_factor) && (x != lbound))
0020    {
0021       BOOST_MATH_INSTRUMENT_VARIABLE(x-1);
0022       return --x;
0023    }
0024    return x;
0025 }
0026
0027 template <class T>
0028 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>&)
0029 {
0030    if((cum < p * fudge_factor) && (x != ubound))
0031    {
0032       BOOST_MATH_INSTRUMENT_VARIABLE(x+1);
0033       return ++x;
0034    }
0035    return x;
0036 }
0037
0038 template <class T>
0039 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>&)
0040 {
0041    if(p >= 0.5)
0042       return round_x_from_p(x, p, cum, fudge_factor, lbound, ubound, policies::discrete_quantile<policies::integer_round_down>());
0043    return round_x_from_p(x, p, cum, fudge_factor, lbound, ubound, policies::discrete_quantile<policies::integer_round_up>());
0044 }
0045
0046 template <class T>
0047 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>&)
0048 {
0049    if(p >= 0.5)
0050       return round_x_from_p(x, p, cum, fudge_factor, lbound, ubound, policies::discrete_quantile<policies::integer_round_up>());
0051    return round_x_from_p(x, p, cum, fudge_factor, lbound, ubound, policies::discrete_quantile<policies::integer_round_down>());
0052 }
0053
0054 template <class T>
0055 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>&)
0056 {
0057    return x;
0058 }
0059
0060 template <class T>
0061 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>&)
0062 {
0063    if((q * fudge_factor > cum) && (x != lbound))
0064    {
0065       BOOST_MATH_INSTRUMENT_VARIABLE(x-1);
0066       return --x;
0067    }
0068    return x;
0069 }
0070
0071 template <class T>
0072 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>&)
0073 {
0074    if((q < cum * fudge_factor) && (x != ubound))
0075    {
0076       BOOST_MATH_INSTRUMENT_VARIABLE(x+1);
0077       return ++x;
0078    }
0079    return x;
0080 }
0081
0082 template <class T>
0083 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>&)
0084 {
0085    if(q < 0.5)
0086       return round_x_from_q(x, q, cum, fudge_factor, lbound, ubound, policies::discrete_quantile<policies::integer_round_down>());
0087    return round_x_from_q(x, q, cum, fudge_factor, lbound, ubound, policies::discrete_quantile<policies::integer_round_up>());
0088 }
0089
0090 template <class T>
0091 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>&)
0092 {
0093    if(q >= 0.5)
0094       return round_x_from_q(x, q, cum, fudge_factor, lbound, ubound, policies::discrete_quantile<policies::integer_round_down>());
0095    return round_x_from_q(x, q, cum, fudge_factor, lbound, ubound, policies::discrete_quantile<policies::integer_round_up>());
0096 }
0097
0098 template <class T>
0099 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>&)
0100 {
0101    return x;
0102 }
0103
0104 template <class T, class Policy>
0105 unsigned hypergeometric_quantile_imp(T p, T q, unsigned r, unsigned n, unsigned N, const Policy& pol)
0106 {
0107 #ifdef _MSC_VER
0108 #  pragma warning(push)
0109 #  pragma warning(disable:4267)
0110 #endif
0111    typedef typename Policy::discrete_quantile_type discrete_quantile_type;
0112    BOOST_MATH_STD_USING
0113    BOOST_FPU_EXCEPTION_GUARD
0114    T result;
0115    T fudge_factor = 1 + tools::epsilon<T>() * ((N <= boost::math::prime(boost::math::max_prime - 1)) ? 50 : 2 * N);
0116    unsigned base = static_cast<unsigned>((std::max)(0, static_cast<int>(n + r) - static_cast<int>(N)));
0117    unsigned lim = (std::min)(r, n);
0118
0119    BOOST_MATH_INSTRUMENT_VARIABLE(p);
0120    BOOST_MATH_INSTRUMENT_VARIABLE(q);
0121    BOOST_MATH_INSTRUMENT_VARIABLE(r);
0122    BOOST_MATH_INSTRUMENT_VARIABLE(n);
0123    BOOST_MATH_INSTRUMENT_VARIABLE(N);
0124    BOOST_MATH_INSTRUMENT_VARIABLE(fudge_factor);
0125    BOOST_MATH_INSTRUMENT_VARIABLE(base);
0126    BOOST_MATH_INSTRUMENT_VARIABLE(lim);
0127
0128    if(p <= 0.5)
0129    {
0130       unsigned x = base;
0131       result = hypergeometric_pdf<T>(x, r, n, N, pol);
0132       T diff = result;
0133       if (diff == 0)
0134       {
0135          ++x;
0136          // We want to skip through x values as fast as we can until we start getting non-zero values,
0137          // otherwise we're just making lots of expensive PDF calls:
0138          T log_pdf = boost::math::lgamma(static_cast<T>(n + 1), pol)
0139             + boost::math::lgamma(static_cast<T>(r + 1), pol)
0140             + boost::math::lgamma(static_cast<T>(N - n + 1), pol)
0141             + boost::math::lgamma(static_cast<T>(N - r + 1), pol)
0142             - boost::math::lgamma(static_cast<T>(N + 1), pol)
0143             - boost::math::lgamma(static_cast<T>(x + 1), pol)
0144             - boost::math::lgamma(static_cast<T>(n - x + 1), pol)
0145             - boost::math::lgamma(static_cast<T>(r - x + 1), pol)
0146             - boost::math::lgamma(static_cast<T>(N - n - r + x + 1), pol);
0147          while (log_pdf < tools::log_min_value<T>())
0148          {
0149             log_pdf += -log(static_cast<T>(x + 1)) + log(static_cast<T>(n - x)) + log(static_cast<T>(r - x)) - log(static_cast<T>(N - n - r + x + 1));
0150             ++x;
0151          }
0152          // By the time we get here, log_pdf may be fairly inaccurate due to
0153          // roundoff errors, get a fresh PDF calculation before proceeding:
0154          diff = hypergeometric_pdf<T>(x, r, n, N, pol);
0155       }
0156       while(result < p)
0157       {
0158          diff = (diff > tools::min_value<T>() * 8)
0159             ? T(n - x) * T(r - x) * diff / (T(x + 1) * T(N + x + 1 - n - r))
0160             : hypergeometric_pdf<T>(x + 1, r, n, N, pol);
0161          if(result + diff / 2 > p)
0162             break;
0163          ++x;
0164          result += diff;
0165 #ifdef BOOST_MATH_INSTRUMENT
0166          if(diff != 0)
0167          {
0168             BOOST_MATH_INSTRUMENT_VARIABLE(x);
0169             BOOST_MATH_INSTRUMENT_VARIABLE(diff);
0170             BOOST_MATH_INSTRUMENT_VARIABLE(result);
0171          }
0172 #endif
0173       }
0174       return round_x_from_p(x, p, result, fudge_factor, base, lim, discrete_quantile_type());
0175    }
0176    else
0177    {
0178       unsigned x = lim;
0179       result = 0;
0180       T diff = hypergeometric_pdf<T>(x, r, n, N, pol);
0181       if (diff == 0)
0182       {
0183          // We want to skip through x values as fast as we can until we start getting non-zero values,
0184          // otherwise we're just making lots of expensive PDF calls:
0185          --x;
0186          T log_pdf = boost::math::lgamma(static_cast<T>(n + 1), pol)
0187             + boost::math::lgamma(static_cast<T>(r + 1), pol)
0188             + boost::math::lgamma(static_cast<T>(N - n + 1), pol)
0189             + boost::math::lgamma(static_cast<T>(N - r + 1), pol)
0190             - boost::math::lgamma(static_cast<T>(N + 1), pol)
0191             - boost::math::lgamma(static_cast<T>(x + 1), pol)
0192             - boost::math::lgamma(static_cast<T>(n - x + 1), pol)
0193             - boost::math::lgamma(static_cast<T>(r - x + 1), pol)
0194             - boost::math::lgamma(static_cast<T>(N - n - r + x + 1), pol);
0195          while (log_pdf < tools::log_min_value<T>())
0196          {
0197             log_pdf += log(static_cast<T>(x)) - log(static_cast<T>(n - x + 1)) - log(static_cast<T>(r - x + 1)) + log(static_cast<T>(N - n - r + x));
0198             --x;
0199          }
0200          // By the time we get here, log_pdf may be fairly inaccurate due to
0201          // roundoff errors, get a fresh PDF calculation before proceeding:
0202          diff = hypergeometric_pdf<T>(x, r, n, N, pol);
0203       }
0204       while(result + diff / 2 < q)
0205       {
0206          result += diff;
0207          diff = (diff > tools::min_value<T>() * 8)
0208             ? x * T(N + x - n - r) * diff / (T(1 + n - x) * T(1 + r - x))
0209             : hypergeometric_pdf<T>(x - 1, r, n, N, pol);
0210          --x;
0211 #ifdef BOOST_MATH_INSTRUMENT
0212          if(diff != 0)
0213          {
0214             BOOST_MATH_INSTRUMENT_VARIABLE(x);
0215             BOOST_MATH_INSTRUMENT_VARIABLE(diff);
0216             BOOST_MATH_INSTRUMENT_VARIABLE(result);
0217          }
0218 #endif
0219       }
0220       return round_x_from_q(x, q, result, fudge_factor, base, lim, discrete_quantile_type());
0221    }
0222 #ifdef _MSC_VER
0223 #  pragma warning(pop)
0224 #endif
0225 }
0226
0227 template <class T, class Policy>
0228 inline unsigned hypergeometric_quantile(T p, T q, unsigned r, unsigned n, unsigned N, const Policy&)
0229 {
0230    BOOST_FPU_EXCEPTION_GUARD
0231    typedef typename tools::promote_args<T>::type result_type;
0232    typedef typename policies::evaluation<result_type, Policy>::type value_type;
0233    typedef typename policies::normalise<
0234       Policy,
0235       policies::promote_float<false>,
0236       policies::promote_double<false>,
0237       policies::assert_undefined<> >::type forwarding_policy;
0238
0239    return detail::hypergeometric_quantile_imp<value_type>(p, q, r, n, N, forwarding_policy());
0240 }
0241
0242 }}} // namespaces
0243
0244 #endif
0245