ceph/src/boost/libs/math/include/boost/math/distributions/detail/hypergeometric_quantile.hpp

   1 // Copyright 2008 John Maddock
   2 //
   3 // Use, modification and distribution are subject to the
   4 // Boost Software License, Version 1.0.
   5 // (See accompanying file LICENSE_1_0.txt
   6 // or copy at http://www.boost.org/LICENSE_1_0.txt)
   7
   8 #ifndef BOOST_MATH_DISTRIBUTIONS_DETAIL_HG_QUANTILE_HPP
   9 #define BOOST_MATH_DISTRIBUTIONS_DETAIL_HG_QUANTILE_HPP
  10
  11 #include <boost/math/policies/error_handling.hpp>
  12 #include <boost/math/distributions/detail/hypergeometric_pdf.hpp>
  13
  14 namespace boost{ namespace math{ namespace detail{
  15
  16 template <class T>
  17 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>&)
  18 {
  19    if((p < cum * fudge_factor) && (x != lbound))
  20    {
  21       BOOST_MATH_INSTRUMENT_VARIABLE(x-1);
  22       return --x;
  23    }
  24    return x;
  25 }
  26
  27 template <class T>
  28 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>&)
  29 {
  30    if((cum < p * fudge_factor) && (x != ubound))
  31    {
  32       BOOST_MATH_INSTRUMENT_VARIABLE(x+1);
  33       return ++x;
  34    }
  35    return x;
  36 }
  37
  38 template <class T>
  39 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>&)
  40 {
  41    if(p >= 0.5)
  42       return round_x_from_p(x, p, cum, fudge_factor, lbound, ubound, policies::discrete_quantile<policies::integer_round_down>());
  43    return round_x_from_p(x, p, cum, fudge_factor, lbound, ubound, policies::discrete_quantile<policies::integer_round_up>());
  44 }
  45
  46 template <class T>
  47 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>&)
  48 {
  49    if(p >= 0.5)
  50       return round_x_from_p(x, p, cum, fudge_factor, lbound, ubound, policies::discrete_quantile<policies::integer_round_up>());
  51    return round_x_from_p(x, p, cum, fudge_factor, lbound, ubound, policies::discrete_quantile<policies::integer_round_down>());
  52 }
  53
  54 template <class T>
  55 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>&)
  56 {
  57    return x;
  58 }
  59
  60 template <class T>
  61 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>&)
  62 {
  63    if((q * fudge_factor > cum) && (x != lbound))
  64    {
  65       BOOST_MATH_INSTRUMENT_VARIABLE(x-1);
  66       return --x;
  67    }
  68    return x;
  69 }
  70
  71 template <class T>
  72 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>&)
  73 {
  74    if((q < cum * fudge_factor) && (x != ubound))
  75    {
  76       BOOST_MATH_INSTRUMENT_VARIABLE(x+1);
  77       return ++x;
  78    }
  79    return x;
  80 }
  81
  82 template <class T>
  83 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>&)
  84 {
  85    if(q < 0.5)
  86       return round_x_from_q(x, q, cum, fudge_factor, lbound, ubound, policies::discrete_quantile<policies::integer_round_down>());
  87    return round_x_from_q(x, q, cum, fudge_factor, lbound, ubound, policies::discrete_quantile<policies::integer_round_up>());
  88 }
  89
  90 template <class T>
  91 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>&)
  92 {
  93    if(q >= 0.5)
  94       return round_x_from_q(x, q, cum, fudge_factor, lbound, ubound, policies::discrete_quantile<policies::integer_round_down>());
  95    return round_x_from_q(x, q, cum, fudge_factor, lbound, ubound, policies::discrete_quantile<policies::integer_round_up>());
  96 }
  97
  98 template <class T>
  99 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>&)
 100 {
 101    return x;
 102 }
 103
 104 template <class T, class Policy>
 105 unsigned hypergeometric_quantile_imp(T p, T q, unsigned r, unsigned n, unsigned N, const Policy& pol)
 106 {
 107 #ifdef BOOST_MSVC
 108 #  pragma warning(push)
 109 #  pragma warning(disable:4267)
 110 #endif
 111    typedef typename Policy::discrete_quantile_type discrete_quantile_type;
 112    BOOST_MATH_STD_USING
 113    BOOST_FPU_EXCEPTION_GUARD
 114    T result;
 115    T fudge_factor = 1 + tools::epsilon<T>() * ((N <= boost::math::prime(boost::math::max_prime - 1)) ? 50 : 2 * N);
 116    unsigned base = static_cast<unsigned>((std::max)(0, (int)(n + r) - (int)(N)));
 117    unsigned lim = (std::min)(r, n);
 118
 119    BOOST_MATH_INSTRUMENT_VARIABLE(p);
 120    BOOST_MATH_INSTRUMENT_VARIABLE(q);
 121    BOOST_MATH_INSTRUMENT_VARIABLE(r);
 122    BOOST_MATH_INSTRUMENT_VARIABLE(n);
 123    BOOST_MATH_INSTRUMENT_VARIABLE(N);
 124    BOOST_MATH_INSTRUMENT_VARIABLE(fudge_factor);
 125    BOOST_MATH_INSTRUMENT_VARIABLE(base);
 126    BOOST_MATH_INSTRUMENT_VARIABLE(lim);
 127
 128    if(p <= 0.5)
 129    {
 130       unsigned x = base;
 131       result = hypergeometric_pdf<T>(x, r, n, N, pol);
 132       T diff = result;
 133       if (diff == 0)
 134       {
 135          ++x;
 136          // We want to skip through x values as fast as we can until we start getting non-zero values,
 137          // otherwise we're just making lots of expensive PDF calls:
 138          T log_pdf = boost::math::lgamma(static_cast<T>(n + 1), pol)
 139             + boost::math::lgamma(static_cast<T>(r + 1), pol)
 140             + boost::math::lgamma(static_cast<T>(N - n + 1), pol)
 141             + boost::math::lgamma(static_cast<T>(N - r + 1), pol)
 142             - boost::math::lgamma(static_cast<T>(N + 1), pol)
 143             - boost::math::lgamma(static_cast<T>(x + 1), pol)
 144             - boost::math::lgamma(static_cast<T>(n - x + 1), pol)
 145             - boost::math::lgamma(static_cast<T>(r - x + 1), pol)
 146             - boost::math::lgamma(static_cast<T>(N - n - r + x + 1), pol);
 147          while (log_pdf < tools::log_min_value<T>())
 148          {
 149             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));
 150             ++x;
 151          }
 152          // By the time we get here, log_pdf may be fairly inaccurate due to
 153          // roundoff errors, get a fresh PDF calculation before proceding:
 154          diff = hypergeometric_pdf<T>(x, r, n, N, pol);
 155       }
 156       while(result < p)
 157       {
 158          diff = (diff > tools::min_value<T>() * 8)
 159             ? T(n - x) * T(r - x) * diff / (T(x + 1) * T(N + x + 1 - n - r))
 160             : hypergeometric_pdf<T>(x + 1, r, n, N, pol);
 161          if(result + diff / 2 > p)
 162             break;
 163          ++x;
 164          result += diff;
 165 #ifdef BOOST_MATH_INSTRUMENT
 166          if(diff != 0)
 167          {
 168             BOOST_MATH_INSTRUMENT_VARIABLE(x);
 169             BOOST_MATH_INSTRUMENT_VARIABLE(diff);
 170             BOOST_MATH_INSTRUMENT_VARIABLE(result);
 171          }
 172 #endif
 173       }
 174       return round_x_from_p(x, p, result, fudge_factor, base, lim, discrete_quantile_type());
 175    }
 176    else
 177    {
 178       unsigned x = lim;
 179       result = 0;
 180       T diff = hypergeometric_pdf<T>(x, r, n, N, pol);
 181       if (diff == 0)
 182       {
 183          // We want to skip through x values as fast as we can until we start getting non-zero values,
 184          // otherwise we're just making lots of expensive PDF calls:
 185          --x;
 186          T log_pdf = boost::math::lgamma(static_cast<T>(n + 1), pol)
 187             + boost::math::lgamma(static_cast<T>(r + 1), pol)
 188             + boost::math::lgamma(static_cast<T>(N - n + 1), pol)
 189             + boost::math::lgamma(static_cast<T>(N - r + 1), pol)
 190             - boost::math::lgamma(static_cast<T>(N + 1), pol)
 191             - boost::math::lgamma(static_cast<T>(x + 1), pol)
 192             - boost::math::lgamma(static_cast<T>(n - x + 1), pol)
 193             - boost::math::lgamma(static_cast<T>(r - x + 1), pol)
 194             - boost::math::lgamma(static_cast<T>(N - n - r + x + 1), pol);
 195          while (log_pdf < tools::log_min_value<T>())
 196          {
 197             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));
 198             --x;
 199          }
 200          // By the time we get here, log_pdf may be fairly inaccurate due to
 201          // roundoff errors, get a fresh PDF calculation before proceding:
 202          diff = hypergeometric_pdf<T>(x, r, n, N, pol);
 203       }
 204       while(result + diff / 2 < q)
 205       {
 206          result += diff;
 207          diff = (diff > tools::min_value<T>() * 8)
 208             ? x * T(N + x - n - r) * diff / (T(1 + n - x) * T(1 + r - x))
 209             : hypergeometric_pdf<T>(x - 1, r, n, N, pol);
 210          --x;
 211 #ifdef BOOST_MATH_INSTRUMENT
 212          if(diff != 0)
 213          {
 214             BOOST_MATH_INSTRUMENT_VARIABLE(x);
 215             BOOST_MATH_INSTRUMENT_VARIABLE(diff);
 216             BOOST_MATH_INSTRUMENT_VARIABLE(result);
 217          }
 218 #endif
 219       }
 220       return round_x_from_q(x, q, result, fudge_factor, base, lim, discrete_quantile_type());
 221    }
 222 #ifdef BOOST_MSVC
 223 #  pragma warning(pop)
 224 #endif
 225 }
 226
 227 template <class T, class Policy>
 228 inline unsigned hypergeometric_quantile(T p, T q, unsigned r, unsigned n, unsigned N, const Policy&)
 229 {
 230    BOOST_FPU_EXCEPTION_GUARD
 231    typedef typename tools::promote_args<T>::type result_type;
 232    typedef typename policies::evaluation<result_type, Policy>::type value_type;
 233    typedef typename policies::normalise<
 234       Policy,
 235       policies::promote_float<false>,
 236       policies::promote_double<false>,
 237       policies::assert_undefined<> >::type forwarding_policy;
 238
 239    return detail::hypergeometric_quantile_imp<value_type>(p, q, r, n, N, forwarding_policy());
 240 }
 241
 242 }}} // namespaces
 243
 244 #endif
 245