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