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());
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());
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    BOOST_MATH_INSTRUMENT_CODE("max_iter = " << max_iter << " count = " << count);
0280    return (r.first + r.second) / 2;
0281 }
0282 //
0283 // Some special routine for rounding up and down:
0284 // We want to check and see if we are very close to an integer, and if so test to see if
0285 // that integer is an exact root of the cdf.  We do this because our root finder only
0286 // guarantees to find *a root*, and there can sometimes be many consecutive floating
0287 // point values which are all roots.  This is especially true if the target probability
0288 // is very close 1.
0289 //
0290 template <class Dist>
0291 inline typename Dist::value_type round_to_floor(const Dist& d, typename Dist::value_type result, typename Dist::value_type p, bool c)
0292 {
0293    BOOST_MATH_STD_USING
0294    typename Dist::value_type cc = ceil(result);
0295    typename Dist::value_type pp = cc <= support(d).second ? c ? cdf(complement(d, cc)) : cdf(d, cc) : 1;
0296    if(pp == p)
0297       result = cc;
0298    else
0299       result = floor(result);
0300    //
0301    // Now find the smallest integer <= result for which we get an exact root:
0302    //
0303    while(result != 0)
0304    {
0305       cc = floor(float_prior(result));
0306       if(cc < support(d).first)
0307          break;
0308       pp = c ? cdf(complement(d, cc)) : cdf(d, cc);
0309       if(c ? pp > p : pp < p)
0310          break;
0311       result = cc;
0312    }
0313
0314    return result;
0315 }
0316
0317 #ifdef _MSC_VER
0318 #pragma warning(push)
0319 #pragma warning(disable:4127)
0320 #endif
0321
0322 template <class Dist>
0323 inline typename Dist::value_type round_to_ceil(const Dist& d, typename Dist::value_type result, typename Dist::value_type p, bool c)
0324 {
0325    BOOST_MATH_STD_USING
0326    typename Dist::value_type cc = floor(result);
0327    typename Dist::value_type pp = cc >= support(d).first ? c ? cdf(complement(d, cc)) : cdf(d, cc) : 0;
0328    if(pp == p)
0329       result = cc;
0330    else
0331       result = ceil(result);
0332    //
0333    // Now find the largest integer >= result for which we get an exact root:
0334    //
0335    while(true)
0336    {
0337       cc = ceil(float_next(result));
0338       if(cc > support(d).second)
0339          break;
0340       pp = c ? cdf(complement(d, cc)) : cdf(d, cc);
0341       if(c ? pp < p : pp > p)
0342          break;
0343       result = cc;
0344    }
0345
0346    return result;
0347 }
0348
0349 #ifdef _MSC_VER
0350 #pragma warning(pop)
0351 #endif
0352 //
0353 // Now finally are the public API functions.
0354 // There is one overload for each policy,
0355 // each one is responsible for selecting the correct
0356 // termination condition, and rounding the result
0357 // to an int where required.
0358 //
0359 template <class Dist>
0360 inline typename Dist::value_type
0361    inverse_discrete_quantile(
0362       const Dist& dist,
0363       typename Dist::value_type p,
0364       bool c,
0365       const typename Dist::value_type& guess,
0366       const typename Dist::value_type& multiplier,
0367       const typename Dist::value_type& adder,
0368       const policies::discrete_quantile<policies::real>&,
0369       std::uintmax_t& max_iter)
0370 {
0371    if(p > 0.5)
0372    {
0373       p = 1 - p;
0374       c = !c;
0375    }
0376    typename Dist::value_type pp = c ? 1 - p : p;
0377    if(pp <= pdf(dist, 0))
0378       return 0;
0379    return do_inverse_discrete_quantile(
0380       dist,
0381       p,
0382       c,
0383       guess,
0384       multiplier,
0385       adder,
0386       tools::eps_tolerance<typename Dist::value_type>(policies::digits<typename Dist::value_type, typename Dist::policy_type>()),
0387       max_iter);
0388 }
0389
0390 template <class Dist>
0391 inline typename Dist::value_type
0392    inverse_discrete_quantile(
0393       const Dist& dist,
0394       const typename Dist::value_type& p,
0395       bool c,
0396       const typename Dist::value_type& guess,
0397       const typename Dist::value_type& multiplier,
0398       const typename Dist::value_type& adder,
0399       const policies::discrete_quantile<policies::integer_round_outwards>&,
0400       std::uintmax_t& max_iter)
0401 {
0402    typedef typename Dist::value_type value_type;
0403    BOOST_MATH_STD_USING
0404    typename Dist::value_type pp = c ? 1 - p : p;
0405    if(pp <= pdf(dist, 0))
0406       return 0;
0407    //
0408    // What happens next depends on whether we're looking for an
0409    // upper or lower quantile:
0410    //
0411    if(pp < 0.5f)
0412       return round_to_floor(dist, do_inverse_discrete_quantile(
0413          dist,
0414          p,
0415          c,
0416          (guess < 1 ? value_type(1) : (value_type)floor(guess)),
0417          multiplier,
0418          adder,
0419          tools::equal_floor(),
0420          max_iter), p, c);
0421    // else:
0422    return round_to_ceil(dist, do_inverse_discrete_quantile(
0423       dist,
0424       p,
0425       c,
0426       (value_type)ceil(guess),
0427       multiplier,
0428       adder,
0429       tools::equal_ceil(),
0430       max_iter), p, c);
0431 }
0432
0433 template <class Dist>
0434 inline typename Dist::value_type
0435    inverse_discrete_quantile(
0436       const Dist& dist,
0437       const typename Dist::value_type& p,
0438       bool c,
0439       const typename Dist::value_type& guess,
0440       const typename Dist::value_type& multiplier,
0441       const typename Dist::value_type& adder,
0442       const policies::discrete_quantile<policies::integer_round_inwards>&,
0443       std::uintmax_t& max_iter)
0444 {
0445    typedef typename Dist::value_type value_type;
0446    BOOST_MATH_STD_USING
0447    typename Dist::value_type pp = c ? 1 - p : p;
0448    if(pp <= pdf(dist, 0))
0449       return 0;
0450    //
0451    // What happens next depends on whether we're looking for an
0452    // upper or lower quantile:
0453    //
0454    if(pp < 0.5f)
0455       return round_to_ceil(dist, do_inverse_discrete_quantile(
0456          dist,
0457          p,
0458          c,
0459          ceil(guess),
0460          multiplier,
0461          adder,
0462          tools::equal_ceil(),
0463          max_iter), p, c);
0464    // else:
0465    return round_to_floor(dist, do_inverse_discrete_quantile(
0466       dist,
0467       p,
0468       c,
0469       (guess < 1 ? value_type(1) : floor(guess)),
0470       multiplier,
0471       adder,
0472       tools::equal_floor(),
0473       max_iter), p, c);
0474 }
0475
0476 template <class Dist>
0477 inline typename Dist::value_type
0478    inverse_discrete_quantile(
0479       const Dist& dist,
0480       const typename Dist::value_type& p,
0481       bool c,
0482       const typename Dist::value_type& guess,
0483       const typename Dist::value_type& multiplier,
0484       const typename Dist::value_type& adder,
0485       const policies::discrete_quantile<policies::integer_round_down>&,
0486       std::uintmax_t& max_iter)
0487 {
0488    typedef typename Dist::value_type value_type;
0489    BOOST_MATH_STD_USING
0490    typename Dist::value_type pp = c ? 1 - p : p;
0491    if(pp <= pdf(dist, 0))
0492       return 0;
0493    return round_to_floor(dist, do_inverse_discrete_quantile(
0494       dist,
0495       p,
0496       c,
0497       (guess < 1 ? value_type(1) : floor(guess)),
0498       multiplier,
0499       adder,
0500       tools::equal_floor(),
0501       max_iter), p, c);
0502 }
0503
0504 template <class Dist>
0505 inline typename Dist::value_type
0506    inverse_discrete_quantile(
0507       const Dist& dist,
0508       const typename Dist::value_type& p,
0509       bool c,
0510       const typename Dist::value_type& guess,
0511       const typename Dist::value_type& multiplier,
0512       const typename Dist::value_type& adder,
0513       const policies::discrete_quantile<policies::integer_round_up>&,
0514       std::uintmax_t& max_iter)
0515 {
0516    BOOST_MATH_STD_USING
0517    typename Dist::value_type pp = c ? 1 - p : p;
0518    if(pp <= pdf(dist, 0))
0519       return 0;
0520    return round_to_ceil(dist, do_inverse_discrete_quantile(
0521       dist,
0522       p,
0523       c,
0524       ceil(guess),
0525       multiplier,
0526       adder,
0527       tools::equal_ceil(),
0528       max_iter), p, c);
0529 }
0530
0531 template <class Dist>
0532 inline typename Dist::value_type
0533    inverse_discrete_quantile(
0534       const Dist& dist,
0535       const typename Dist::value_type& p,
0536       bool c,
0537       const typename Dist::value_type& guess,
0538       const typename Dist::value_type& multiplier,
0539       const typename Dist::value_type& adder,
0540       const policies::discrete_quantile<policies::integer_round_nearest>&,
0541       std::uintmax_t& max_iter)
0542 {
0543    typedef typename Dist::value_type value_type;
0544    BOOST_MATH_STD_USING
0545    typename Dist::value_type pp = c ? 1 - p : p;
0546    if(pp <= pdf(dist, 0))
0547       return 0;
0548    //
0549    // Note that we adjust the guess to the nearest half-integer:
0550    // this increase the chances that we will bracket the root
0551    // with two results that both round to the same integer quickly.
0552    //
0553    return round_to_floor(dist, do_inverse_discrete_quantile(
0554       dist,
0555       p,
0556       c,
0557       (guess < 0.5f ? value_type(1.5f) : floor(guess + 0.5f) + 0.5f),
0558       multiplier,
0559       adder,
0560       tools::equal_nearest_integer(),
0561       max_iter) + 0.5f, p, c);
0562 }
0563
0564 }}} // namespaces
0565
0566 #endif // BOOST_MATH_DISTRIBUTIONS_DETAIL_INV_DISCRETE_QUANTILE
0567