distributions/detail/inv_discrete_quantile.hpp

0001 //  Copyright John Maddock 2007.
0002 //  Use, modification and distribution are subject to the
0003 //  Boost Software License, Version 1.0. (See accompanying file
0004 //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
0005
0006 #ifndef BOOST_MATH_DISTRIBUTIONS_DETAIL_INV_DISCRETE_QUANTILE
0007 #define BOOST_MATH_DISTRIBUTIONS_DETAIL_INV_DISCRETE_QUANTILE
0008
0009 #include <boost/math/tools/config.hpp>
0010 #include <boost/math/tools/cstdint.hpp>
0011 #include <boost/math/tools/precision.hpp>
0012 #include <boost/math/tools/toms748_solve.hpp>
0013 #include <boost/math/tools/tuple.hpp>
0014
0015 namespace boost{ namespace math{ namespace detail{
0016
0017 //
0018 // Functor for root finding algorithm:
0019 //
0020 template <class Dist>
0021 struct distribution_quantile_finder
0022 {
0023    typedef typename Dist::value_type value_type;
0024    typedef typename Dist::policy_type policy_type;
0025
0026    BOOST_MATH_GPU_ENABLED distribution_quantile_finder(const Dist d, value_type p, bool c)
0027       : dist(d), target(p), comp(c) {}
0028
0029    BOOST_MATH_GPU_ENABLED value_type operator()(value_type const& x)
0030    {
0031       return comp ? value_type(target - cdf(complement(dist, x))) : value_type(cdf(dist, x) - target);
0032    }
0033
0034 private:
0035    Dist dist;
0036    value_type target;
0037    bool comp;
0038 };
0039 //
0040 // The purpose of adjust_bounds, is to toggle the last bit of the
0041 // range so that both ends round to the same integer, if possible.
0042 // If they do both round the same then we terminate the search
0043 // for the root *very* quickly when finding an integer result.
0044 // At the point that this function is called we know that "a" is
0045 // below the root and "b" above it, so this change can not result
0046 // in the root no longer being bracketed.
0047 //
0048 template <class Real, class Tol>
0049 BOOST_MATH_GPU_ENABLED void adjust_bounds(Real& /* a */, Real& /* b */, Tol const& /* tol */){}
0050
0051 template <class Real>
0052 BOOST_MATH_GPU_ENABLED void adjust_bounds(Real& /* a */, Real& b, tools::equal_floor const& /* tol */)
0053 {
0054    BOOST_MATH_STD_USING
0055    b -= tools::epsilon<Real>() * b;
0056 }
0057
0058 template <class Real>
0059 BOOST_MATH_GPU_ENABLED void adjust_bounds(Real& a, Real& /* b */, tools::equal_ceil const& /* tol */)
0060 {
0061    BOOST_MATH_STD_USING
0062    a += tools::epsilon<Real>() * a;
0063 }
0064
0065 template <class Real>
0066 BOOST_MATH_GPU_ENABLED void adjust_bounds(Real& a, Real& b, tools::equal_nearest_integer const& /* tol */)
0067 {
0068    BOOST_MATH_STD_USING
0069    a += tools::epsilon<Real>() * a;
0070    b -= tools::epsilon<Real>() * b;
0071 }
0072 //
0073 // This is where all the work is done:
0074 //
0075 template <class Dist, class Tolerance>
0076 BOOST_MATH_GPU_ENABLED typename Dist::value_type
0077    do_inverse_discrete_quantile(
0078       const Dist& dist,
0079       const typename Dist::value_type& p,
0080       bool comp,
0081       typename Dist::value_type guess,
0082       const typename Dist::value_type& multiplier,
0083       typename Dist::value_type adder,
0084       const Tolerance& tol,
0085       boost::math::uintmax_t& max_iter)
0086 {
0087    typedef typename Dist::value_type value_type;
0088    typedef typename Dist::policy_type policy_type;
0089
0090    constexpr auto function = "boost::math::do_inverse_discrete_quantile<%1%>";
0091
0092    BOOST_MATH_STD_USING
0093
0094    distribution_quantile_finder<Dist> f(dist, p, comp);
0095    //
0096    // Max bounds of the distribution:
0097    //
0098    value_type min_bound, max_bound;
0099    boost::math::tie(min_bound, max_bound) = support(dist);
0100
0101    if(guess > max_bound)
0102       guess = max_bound;
0103    if(guess < min_bound)
0104       guess = min_bound;
0105
0106    value_type fa = f(guess);
0107    boost::math::uintmax_t count = max_iter - 1;
0108    value_type fb(fa), a(guess), b =0; // Compiler warning C4701: potentially uninitialized local variable 'b' used
0109
0110    if(fa == 0)
0111       return guess;
0112
0113    //
0114    // For small expected results, just use a linear search:
0115    //
0116    if(guess < 10)
0117    {
0118       b = a;
0119       while((a < 10) && (fa * fb >= 0))
0120       {
0121          if(fb <= 0)
0122          {
0123             a = b;
0124             b = a + 1;
0125             if(b > max_bound)
0126                b = max_bound;
0127             fb = f(b);
0128             --count;
0129             if(fb == 0)
0130                return b;
0131             if(a == b)
0132                return b; // can't go any higher!
0133          }
0134          else
0135          {
0136             b = a;
0137             a = BOOST_MATH_GPU_SAFE_MAX(value_type(b - 1), value_type(0));
0138             if(a < min_bound)
0139                a = min_bound;
0140             fa = f(a);
0141             --count;
0142             if(fa == 0)
0143                return a;
0144             if(a == b)
0145                return a;  //  We can't go any lower than this!
0146          }
0147       }
0148    }
0149    //
0150    // Try and bracket using a couple of additions first,
0151    // we're assuming that "guess" is likely to be accurate
0152    // to the nearest int or so:
0153    //
0154    else if((adder != 0) && (a + adder != a))
0155    {
0156       //
0157       // If we're looking for a large result, then bump "adder" up
0158       // by a bit to increase our chances of bracketing the root:
0159       //
0160       //adder = BOOST_MATH_GPU_SAFE_MAX(adder, 0.001f * guess);
0161       if(fa < 0)
0162       {
0163          b = a + adder;
0164          if(b > max_bound)
0165             b = max_bound;
0166       }
0167       else
0168       {
0169          b = BOOST_MATH_GPU_SAFE_MAX(value_type(a - adder), value_type(0));
0170          if(b < min_bound)
0171             b = min_bound;
0172       }
0173       fb = f(b);
0174       --count;
0175       if(fb == 0)
0176          return b;
0177       if(count && (fa * fb >= 0))
0178       {
0179          //
0180          // We didn't bracket the root, try
0181          // once more:
0182          //
0183          a = b;
0184          fa = fb;
0185          if(fa < 0)
0186          {
0187             b = a + adder;
0188             if(b > max_bound)
0189                b = max_bound;
0190          }
0191          else
0192          {
0193             b = BOOST_MATH_GPU_SAFE_MAX(value_type(a - adder), value_type(0));
0194             if(b < min_bound)
0195                b = min_bound;
0196          }
0197          fb = f(b);
0198          --count;
0199       }
0200       if(a > b)
0201       {
0202          BOOST_MATH_GPU_SAFE_SWAP(a, b);
0203          BOOST_MATH_GPU_SAFE_SWAP(fa, fb);
0204       }
0205    }
0206    //
0207    // If the root hasn't been bracketed yet, try again
0208    // using the multiplier this time:
0209    //
0210    if((boost::math::sign)(fb) == (boost::math::sign)(fa))
0211    {
0212       if(fa < 0)
0213       {
0214          //
0215          // Zero is to the right of x2, so walk upwards
0216          // until we find it:
0217          //
0218          while(((boost::math::sign)(fb) == (boost::math::sign)(fa)) && (a != b))
0219          {
0220             if(count == 0)
0221                return policies::raise_evaluation_error(function, "Unable to bracket root, last nearest value was %1%", b, policy_type()); // LCOV_EXCL_LINE
0222             a = b;
0223             fa = fb;
0224             b *= multiplier;
0225             if(b > max_bound)
0226                b = max_bound;
0227             fb = f(b);
0228             --count;
0229             BOOST_MATH_INSTRUMENT_CODE("a = " << a << " b = " << b << " fa = " << fa << " fb = " << fb << " count = " << count);
0230          }
0231       }
0232       else
0233       {
0234          //
0235          // Zero is to the left of a, so walk downwards
0236          // until we find it:
0237          //
0238          while(((boost::math::sign)(fb) == (boost::math::sign)(fa)) && (a != b))
0239          {
0240             if(fabs(a) < tools::min_value<value_type>())
0241             {
0242                // Escape route just in case the answer is zero!
0243                max_iter -= count;
0244                max_iter += 1;
0245                return 0;
0246             }
0247             if(count == 0)
0248                return policies::raise_evaluation_error(function, "Unable to bracket root, last nearest value was %1%", a, policy_type()); // LCOV_EXCL_LINE
0249             b = a;
0250             fb = fa;
0251             a /= multiplier;
0252             if(a < min_bound)
0253                a = min_bound;
0254             fa = f(a);
0255             --count;
0256             BOOST_MATH_INSTRUMENT_CODE("a = " << a << " b = " << b << " fa = " << fa << " fb = " << fb << " count = " << count);
0257          }
0258       }
0259    }
0260    max_iter -= count;
0261    if(fa == 0)
0262       return a;
0263    if(fb == 0)
0264       return b;
0265    if(a == b)
0266       return b;  // Ran out of bounds trying to bracket - there is no answer!
0267    //
0268    // Adjust bounds so that if we're looking for an integer
0269    // result, then both ends round the same way:
0270    //
0271    adjust_bounds(a, b, tol);
0272    //
0273    // We don't want zero or denorm lower bounds:
0274    //
0275    if(a < tools::min_value<value_type>())
0276       a = tools::min_value<value_type>();
0277    //
0278    // Go ahead and find the root:
0279    //
0280    boost::math::pair<value_type, value_type> r = toms748_solve(f, a, b, fa, fb, tol, count, policy_type());
0281    max_iter += count;
0282    if (max_iter >= policies::get_max_root_iterations<policy_type>())
0283    {
0284       return policies::raise_evaluation_error<value_type>(function, "Unable to locate solution in a reasonable time:" // LCOV_EXCL_LINE
0285          " either there is no answer to quantile or the answer is infinite.  Current best guess is %1%", r.first, policy_type()); // LCOV_EXCL_LINE
0286    }
0287    BOOST_MATH_INSTRUMENT_CODE("max_iter = " << max_iter << " count = " << count);
0288    return (r.first + r.second) / 2;
0289 }
0290 //
0291 // Some special routine for rounding up and down:
0292 // We want to check and see if we are very close to an integer, and if so test to see if
0293 // that integer is an exact root of the cdf.  We do this because our root finder only
0294 // guarantees to find *a root*, and there can sometimes be many consecutive floating
0295 // point values which are all roots.  This is especially true if the target probability
0296 // is very close 1.
0297 //
0298 template <class Dist>
0299 BOOST_MATH_GPU_ENABLED inline typename Dist::value_type round_to_floor(const Dist& d, typename Dist::value_type result, typename Dist::value_type p, bool c)
0300 {
0301    BOOST_MATH_STD_USING
0302    typename Dist::value_type cc = ceil(result);
0303    typename Dist::value_type pp = cc <= support(d).second ? c ? cdf(complement(d, cc)) : cdf(d, cc) : 1;
0304    if(pp == p)
0305       result = cc;
0306    else
0307       result = floor(result);
0308    //
0309    // Now find the smallest integer <= result for which we get an exact root:
0310    //
0311    while(result != 0)
0312    {
0313       #ifdef BOOST_MATH_HAS_GPU_SUPPORT
0314       cc = floor(::nextafter(result, -tools::max_value<typename Dist::value_type>()));
0315       #else
0316       cc = floor(float_prior(result));
0317       #endif
0318       if(cc < support(d).first)
0319          break;
0320       pp = c ? cdf(complement(d, cc)) : cdf(d, cc);
0321       if(c ? pp > p : pp < p)
0322          break;
0323       result = cc;
0324    }
0325
0326    return result;
0327 }
0328
0329 #ifdef _MSC_VER
0330 #pragma warning(push)
0331 #pragma warning(disable:4127)
0332 #endif
0333
0334 template <class Dist>
0335 BOOST_MATH_GPU_ENABLED inline typename Dist::value_type round_to_ceil(const Dist& d, typename Dist::value_type result, typename Dist::value_type p, bool c)
0336 {
0337    BOOST_MATH_STD_USING
0338    typename Dist::value_type cc = floor(result);
0339    typename Dist::value_type pp = cc >= support(d).first ? c ? cdf(complement(d, cc)) : cdf(d, cc) : 0;
0340    if(pp == p)
0341       result = cc;
0342    else
0343       result = ceil(result);
0344    //
0345    // Now find the largest integer >= result for which we get an exact root:
0346    //
0347    while(true)
0348    {
0349       #ifdef BOOST_MATH_HAS_GPU_SUPPORT
0350       cc = ceil(::nextafter(result, tools::max_value<typename Dist::value_type>()));
0351       #else
0352       cc = ceil(float_next(result));
0353       #endif
0354       if(cc > support(d).second)
0355          break;
0356       pp = c ? cdf(complement(d, cc)) : cdf(d, cc);
0357       if(c ? pp < p : pp > p)
0358          break;
0359       result = cc;
0360    }
0361
0362    return result;
0363 }
0364
0365 #ifdef _MSC_VER
0366 #pragma warning(pop)
0367 #endif
0368 //
0369 // Now finally are the public API functions.
0370 // There is one overload for each policy,
0371 // each one is responsible for selecting the correct
0372 // termination condition, and rounding the result
0373 // to an int where required.
0374 //
0375 template <class Dist>
0376 BOOST_MATH_GPU_ENABLED inline typename Dist::value_type
0377    inverse_discrete_quantile(
0378       const Dist& dist,
0379       typename Dist::value_type p,
0380       bool c,
0381       const typename Dist::value_type& guess,
0382       const typename Dist::value_type& multiplier,
0383       const typename Dist::value_type& adder,
0384       const policies::discrete_quantile<policies::real>&,
0385       boost::math::uintmax_t& max_iter)
0386 {
0387    if(p > 0.5)
0388    {
0389       p = 1 - p;
0390       c = !c;
0391    }
0392    typename Dist::value_type pp = c ? 1 - p : p;
0393    if(pp <= pdf(dist, 0))
0394       return 0;
0395    return do_inverse_discrete_quantile(
0396       dist,
0397       p,
0398       c,
0399       guess,
0400       multiplier,
0401       adder,
0402       tools::eps_tolerance<typename Dist::value_type>(policies::digits<typename Dist::value_type, typename Dist::policy_type>()),
0403       max_iter);
0404 }
0405
0406 template <class Dist>
0407 BOOST_MATH_GPU_ENABLED inline typename Dist::value_type
0408    inverse_discrete_quantile(
0409       const Dist& dist,
0410       const typename Dist::value_type& p,
0411       bool c,
0412       const typename Dist::value_type& guess,
0413       const typename Dist::value_type& multiplier,
0414       const typename Dist::value_type& adder,
0415       const policies::discrete_quantile<policies::integer_round_outwards>&,
0416       boost::math::uintmax_t& max_iter)
0417 {
0418    typedef typename Dist::value_type value_type;
0419    BOOST_MATH_STD_USING
0420    typename Dist::value_type pp = c ? 1 - p : p;
0421    if(pp <= pdf(dist, 0))
0422       return 0;
0423    //
0424    // What happens next depends on whether we're looking for an
0425    // upper or lower quantile:
0426    //
0427    if(pp < 0.5f)
0428       return round_to_floor(dist, do_inverse_discrete_quantile(
0429          dist,
0430          p,
0431          c,
0432          (guess < 1 ? value_type(1) : (value_type)floor(guess)),
0433          multiplier,
0434          adder,
0435          tools::equal_floor(),
0436          max_iter), p, c);
0437    // else:
0438    return round_to_ceil(dist, do_inverse_discrete_quantile(
0439       dist,
0440       p,
0441       c,
0442       (value_type)ceil(guess),
0443       multiplier,
0444       adder,
0445       tools::equal_ceil(),
0446       max_iter), p, c);
0447 }
0448
0449 template <class Dist>
0450 BOOST_MATH_GPU_ENABLED inline typename Dist::value_type
0451    inverse_discrete_quantile(
0452       const Dist& dist,
0453       const typename Dist::value_type& p,
0454       bool c,
0455       const typename Dist::value_type& guess,
0456       const typename Dist::value_type& multiplier,
0457       const typename Dist::value_type& adder,
0458       const policies::discrete_quantile<policies::integer_round_inwards>&,
0459       boost::math::uintmax_t& max_iter)
0460 {
0461    typedef typename Dist::value_type value_type;
0462    BOOST_MATH_STD_USING
0463    typename Dist::value_type pp = c ? 1 - p : p;
0464    if(pp <= pdf(dist, 0))
0465       return 0;
0466    //
0467    // What happens next depends on whether we're looking for an
0468    // upper or lower quantile:
0469    //
0470    if(pp < 0.5f)
0471       return round_to_ceil(dist, do_inverse_discrete_quantile(
0472          dist,
0473          p,
0474          c,
0475          ceil(guess),
0476          multiplier,
0477          adder,
0478          tools::equal_ceil(),
0479          max_iter), p, c);
0480    // else:
0481    return round_to_floor(dist, do_inverse_discrete_quantile(
0482       dist,
0483       p,
0484       c,
0485       (guess < 1 ? value_type(1) : floor(guess)),
0486       multiplier,
0487       adder,
0488       tools::equal_floor(),
0489       max_iter), p, c);
0490 }
0491
0492 template <class Dist>
0493 BOOST_MATH_GPU_ENABLED inline typename Dist::value_type
0494    inverse_discrete_quantile(
0495       const Dist& dist,
0496       const typename Dist::value_type& p,
0497       bool c,
0498       const typename Dist::value_type& guess,
0499       const typename Dist::value_type& multiplier,
0500       const typename Dist::value_type& adder,
0501       const policies::discrete_quantile<policies::integer_round_down>&,
0502       boost::math::uintmax_t& max_iter)
0503 {
0504    typedef typename Dist::value_type value_type;
0505    BOOST_MATH_STD_USING
0506    typename Dist::value_type pp = c ? 1 - p : p;
0507    if(pp <= pdf(dist, 0))
0508       return 0;
0509    return round_to_floor(dist, do_inverse_discrete_quantile(
0510       dist,
0511       p,
0512       c,
0513       (guess < 1 ? value_type(1) : floor(guess)),
0514       multiplier,
0515       adder,
0516       tools::equal_floor(),
0517       max_iter), p, c);
0518 }
0519
0520 template <class Dist>
0521 BOOST_MATH_GPU_ENABLED inline typename Dist::value_type
0522    inverse_discrete_quantile(
0523       const Dist& dist,
0524       const typename Dist::value_type& p,
0525       bool c,
0526       const typename Dist::value_type& guess,
0527       const typename Dist::value_type& multiplier,
0528       const typename Dist::value_type& adder,
0529       const policies::discrete_quantile<policies::integer_round_up>&,
0530       boost::math::uintmax_t& max_iter)
0531 {
0532    BOOST_MATH_STD_USING
0533    typename Dist::value_type pp = c ? 1 - p : p;
0534    if(pp <= pdf(dist, 0))
0535       return 0;
0536    return round_to_ceil(dist, do_inverse_discrete_quantile(
0537       dist,
0538       p,
0539       c,
0540       ceil(guess),
0541       multiplier,
0542       adder,
0543       tools::equal_ceil(),
0544       max_iter), p, c);
0545 }
0546
0547 template <class Dist>
0548 BOOST_MATH_GPU_ENABLED inline typename Dist::value_type
0549    inverse_discrete_quantile(
0550       const Dist& dist,
0551       const typename Dist::value_type& p,
0552       bool c,
0553       const typename Dist::value_type& guess,
0554       const typename Dist::value_type& multiplier,
0555       const typename Dist::value_type& adder,
0556       const policies::discrete_quantile<policies::integer_round_nearest>&,
0557       boost::math::uintmax_t& max_iter)
0558 {
0559    typedef typename Dist::value_type value_type;
0560    BOOST_MATH_STD_USING
0561    typename Dist::value_type pp = c ? 1 - p : p;
0562    if(pp <= pdf(dist, 0))
0563       return 0;
0564    //
0565    // Note that we adjust the guess to the nearest half-integer:
0566    // this increase the chances that we will bracket the root
0567    // with two results that both round to the same integer quickly.
0568    //
0569    return round_to_floor(dist, do_inverse_discrete_quantile(
0570       dist,
0571       p,
0572       c,
0573       (guess < 0.5f ? value_type(1.5f) : floor(guess + 0.5f) + 0.5f),
0574       multiplier,
0575       adder,
0576       tools::equal_nearest_integer(),
0577       max_iter) + 0.5f, p, c);
0578 }
0579
0580 }}} // namespaces
0581
0582 #endif // BOOST_MATH_DISTRIBUTIONS_DETAIL_INV_DISCRETE_QUANTILE
0583