ceph/src/boost/boost/math/distributions/detail/inv_discrete_quantile.hpp

   1 //  Copyright John Maddock 2007.
   2 //  Use, modification and distribution are subject to the
   3 //  Boost Software License, Version 1.0. (See accompanying file
   4 //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
   5
   6 #ifndef BOOST_MATH_DISTRIBUTIONS_DETAIL_INV_DISCRETE_QUANTILE
   7 #define BOOST_MATH_DISTRIBUTIONS_DETAIL_INV_DISCRETE_QUANTILE
   8
   9 #include <algorithm>
  10
  11 namespace boost{ namespace math{ namespace detail{
  12
  13 //
  14 // Functor for root finding algorithm:
  15 //
  16 template <class Dist>
  17 struct distribution_quantile_finder
  18 {
  19    typedef typename Dist::value_type value_type;
  20    typedef typename Dist::policy_type policy_type;
  21
  22    distribution_quantile_finder(const Dist d, value_type p, bool c)
  23       : dist(d), target(p), comp(c) {}
  24
  25    value_type operator()(value_type const& x)
  26    {
  27       return comp ? value_type(target - cdf(complement(dist, x))) : value_type(cdf(dist, x) - target);
  28    }
  29
  30 private:
  31    Dist dist;
  32    value_type target;
  33    bool comp;
  34 };
  35 //
  36 // The purpose of adjust_bounds, is to toggle the last bit of the
  37 // range so that both ends round to the same integer, if possible.
  38 // If they do both round the same then we terminate the search
  39 // for the root *very* quickly when finding an integer result.
  40 // At the point that this function is called we know that "a" is
  41 // below the root and "b" above it, so this change can not result
  42 // in the root no longer being bracketed.
  43 //
  44 template <class Real, class Tol>
  45 void adjust_bounds(Real& /* a */, Real& /* b */, Tol const& /* tol */){}
  46
  47 template <class Real>
  48 void adjust_bounds(Real& /* a */, Real& b, tools::equal_floor const& /* tol */)
  49 {
  50    BOOST_MATH_STD_USING
  51    b -= tools::epsilon<Real>() * b;
  52 }
  53
  54 template <class Real>
  55 void adjust_bounds(Real& a, Real& /* b */, tools::equal_ceil const& /* tol */)
  56 {
  57    BOOST_MATH_STD_USING
  58    a += tools::epsilon<Real>() * a;
  59 }
  60
  61 template <class Real>
  62 void adjust_bounds(Real& a, Real& b, tools::equal_nearest_integer const& /* tol */)
  63 {
  64    BOOST_MATH_STD_USING
  65    a += tools::epsilon<Real>() * a;
  66    b -= tools::epsilon<Real>() * b;
  67 }
  68 //
  69 // This is where all the work is done:
  70 //
  71 template <class Dist, class Tolerance>
  72 typename Dist::value_type
  73    do_inverse_discrete_quantile(
  74       const Dist& dist,
  75       const typename Dist::value_type& p,
  76       bool comp,
  77       typename Dist::value_type guess,
  78       const typename Dist::value_type& multiplier,
  79       typename Dist::value_type adder,
  80       const Tolerance& tol,
  81       boost::uintmax_t& max_iter)
  82 {
  83    typedef typename Dist::value_type value_type;
  84    typedef typename Dist::policy_type policy_type;
  85
  86    static const char* function = "boost::math::do_inverse_discrete_quantile<%1%>";
  87
  88    BOOST_MATH_STD_USING
  89
  90    distribution_quantile_finder<Dist> f(dist, p, comp);
  91    //
  92    // Max bounds of the distribution:
  93    //
  94    value_type min_bound, max_bound;
  95    boost::math::tie(min_bound, max_bound) = support(dist);
  96
  97    if(guess > max_bound)
  98       guess = max_bound;
  99    if(guess < min_bound)
 100       guess = min_bound;
 101
 102    value_type fa = f(guess);
 103    boost::uintmax_t count = max_iter - 1;
 104    value_type fb(fa), a(guess), b =0; // Compiler warning C4701: potentially uninitialized local variable 'b' used
 105
 106    if(fa == 0)
 107       return guess;
 108
 109    //
 110    // For small expected results, just use a linear search:
 111    //
 112    if(guess < 10)
 113    {
 114       b = a;
 115       while((a < 10) && (fa * fb >= 0))
 116       {
 117          if(fb <= 0)
 118          {
 119             a = b;
 120             b = a + 1;
 121             if(b > max_bound)
 122                b = max_bound;
 123             fb = f(b);
 124             --count;
 125             if(fb == 0)
 126                return b;
 127             if(a == b)
 128                return b; // can't go any higher!
 129          }
 130          else
 131          {
 132             b = a;
 133             a = (std::max)(value_type(b - 1), value_type(0));
 134             if(a < min_bound)
 135                a = min_bound;
 136             fa = f(a);
 137             --count;
 138             if(fa == 0)
 139                return a;
 140             if(a == b)
 141                return a;  //  We can't go any lower than this!
 142          }
 143       }
 144    }
 145    //
 146    // Try and bracket using a couple of additions first,
 147    // we're assuming that "guess" is likely to be accurate
 148    // to the nearest int or so:
 149    //
 150    else if(adder != 0)
 151    {
 152       //
 153       // If we're looking for a large result, then bump "adder" up
 154       // by a bit to increase our chances of bracketing the root:
 155       //
 156       //adder = (std::max)(adder, 0.001f * guess);
 157       if(fa < 0)
 158       {
 159          b = a + adder;
 160          if(b > max_bound)
 161             b = max_bound;
 162       }
 163       else
 164       {
 165          b = (std::max)(value_type(a - adder), value_type(0));
 166          if(b < min_bound)
 167             b = min_bound;
 168       }
 169       fb = f(b);
 170       --count;
 171       if(fb == 0)
 172          return b;
 173       if(count && (fa * fb >= 0))
 174       {
 175          //
 176          // We didn't bracket the root, try
 177          // once more:
 178          //
 179          a = b;
 180          fa = fb;
 181          if(fa < 0)
 182          {
 183             b = a + adder;
 184             if(b > max_bound)
 185                b = max_bound;
 186          }
 187          else
 188          {
 189             b = (std::max)(value_type(a - adder), value_type(0));
 190             if(b < min_bound)
 191                b = min_bound;
 192          }
 193          fb = f(b);
 194          --count;
 195       }
 196       if(a > b)
 197       {
 198          using std::swap;
 199          swap(a, b);
 200          swap(fa, fb);
 201       }
 202    }
 203    //
 204    // If the root hasn't been bracketed yet, try again
 205    // using the multiplier this time:
 206    //
 207    if((boost::math::sign)(fb) == (boost::math::sign)(fa))
 208    {
 209       if(fa < 0)
 210       {
 211          //
 212          // Zero is to the right of x2, so walk upwards
 213          // until we find it:
 214          //
 215          while(((boost::math::sign)(fb) == (boost::math::sign)(fa)) && (a != b))
 216          {
 217             if(count == 0)
 218                return policies::raise_evaluation_error(function, "Unable to bracket root, last nearest value was %1%", b, policy_type());
 219             a = b;
 220             fa = fb;
 221             b *= multiplier;
 222             if(b > max_bound)
 223                b = max_bound;
 224             fb = f(b);
 225             --count;
 226             BOOST_MATH_INSTRUMENT_CODE("a = " << a << " b = " << b << " fa = " << fa << " fb = " << fb << " count = " << count);
 227          }
 228       }
 229       else
 230       {
 231          //
 232          // Zero is to the left of a, so walk downwards
 233          // until we find it:
 234          //
 235          while(((boost::math::sign)(fb) == (boost::math::sign)(fa)) && (a != b))
 236          {
 237             if(fabs(a) < tools::min_value<value_type>())
 238             {
 239                // Escape route just in case the answer is zero!
 240                max_iter -= count;
 241                max_iter += 1;
 242                return 0;
 243             }
 244             if(count == 0)
 245                return policies::raise_evaluation_error(function, "Unable to bracket root, last nearest value was %1%", a, policy_type());
 246             b = a;
 247             fb = fa;
 248             a /= multiplier;
 249             if(a < min_bound)
 250                a = min_bound;
 251             fa = f(a);
 252             --count;
 253             BOOST_MATH_INSTRUMENT_CODE("a = " << a << " b = " << b << " fa = " << fa << " fb = " << fb << " count = " << count);
 254          }
 255       }
 256    }
 257    max_iter -= count;
 258    if(fa == 0)
 259       return a;
 260    if(fb == 0)
 261       return b;
 262    if(a == b)
 263       return b;  // Ran out of bounds trying to bracket - there is no answer!
 264    //
 265    // Adjust bounds so that if we're looking for an integer
 266    // result, then both ends round the same way:
 267    //
 268    adjust_bounds(a, b, tol);
 269    //
 270    // We don't want zero or denorm lower bounds:
 271    //
 272    if(a < tools::min_value<value_type>())
 273       a = tools::min_value<value_type>();
 274    //
 275    // Go ahead and find the root:
 276    //
 277    std::pair<value_type, value_type> r = toms748_solve(f, a, b, fa, fb, tol, count, policy_type());
 278    max_iter += count;
 279    BOOST_MATH_INSTRUMENT_CODE("max_iter = " << max_iter << " count = " << count);
 280    return (r.first + r.second) / 2;
 281 }
 282 //
 283 // Some special routine for rounding up and down:
 284 // We want to check and see if we are very close to an integer, and if so test to see if
 285 // that integer is an exact root of the cdf.  We do this because our root finder only
 286 // guarantees to find *a root*, and there can sometimes be many consecutive floating
 287 // point values which are all roots.  This is especially true if the target probability
 288 // is very close 1.
 289 //
 290 template <class Dist>
 291 inline typename Dist::value_type round_to_floor(const Dist& d, typename Dist::value_type result, typename Dist::value_type p, bool c)
 292 {
 293    BOOST_MATH_STD_USING
 294    typename Dist::value_type cc = ceil(result);
 295    typename Dist::value_type pp = cc <= support(d).second ? c ? cdf(complement(d, cc)) : cdf(d, cc) : 1;
 296    if(pp == p)
 297       result = cc;
 298    else
 299       result = floor(result);
 300    //
 301    // Now find the smallest integer <= result for which we get an exact root:
 302    //
 303    while(result != 0)
 304    {
 305       cc = result - 1;
 306       if(cc < support(d).first)
 307          break;
 308       pp = c ? cdf(complement(d, cc)) : cdf(d, cc);
 309       if(pp == p)
 310          result = cc;
 311       else if(c ? pp > p : pp < p)
 312          break;
 313       result -= 1;
 314    }
 315
 316    return result;
 317 }
 318
 319 #ifdef BOOST_MSVC
 320 #pragma warning(push)
 321 #pragma warning(disable:4127)
 322 #endif
 323
 324 template <class Dist>
 325 inline typename Dist::value_type round_to_ceil(const Dist& d, typename Dist::value_type result, typename Dist::value_type p, bool c)
 326 {
 327    BOOST_MATH_STD_USING
 328    typename Dist::value_type cc = floor(result);
 329    typename Dist::value_type pp = cc >= support(d).first ? c ? cdf(complement(d, cc)) : cdf(d, cc) : 0;
 330    if(pp == p)
 331       result = cc;
 332    else
 333       result = ceil(result);
 334    //
 335    // Now find the largest integer >= result for which we get an exact root:
 336    //
 337    while(true)
 338    {
 339       cc = result + 1;
 340       if(cc > support(d).second)
 341          break;
 342       pp = c ? cdf(complement(d, cc)) : cdf(d, cc);
 343       if(pp == p)
 344          result = cc;
 345       else if(c ? pp < p : pp > p)
 346          break;
 347       result += 1;
 348    }
 349
 350    return result;
 351 }
 352
 353 #ifdef BOOST_MSVC
 354 #pragma warning(pop)
 355 #endif
 356 //
 357 // Now finally are the public API functions.
 358 // There is one overload for each policy,
 359 // each one is responsible for selecting the correct
 360 // termination condition, and rounding the result
 361 // to an int where required.
 362 //
 363 template <class Dist>
 364 inline typename Dist::value_type
 365    inverse_discrete_quantile(
 366       const Dist& dist,
 367       typename Dist::value_type p,
 368       bool c,
 369       const typename Dist::value_type& guess,
 370       const typename Dist::value_type& multiplier,
 371       const typename Dist::value_type& adder,
 372       const policies::discrete_quantile<policies::real>&,
 373       boost::uintmax_t& max_iter)
 374 {
 375    if(p > 0.5)
 376    {
 377       p = 1 - p;
 378       c = !c;
 379    }
 380    typename Dist::value_type pp = c ? 1 - p : p;
 381    if(pp <= pdf(dist, 0))
 382       return 0;
 383    return do_inverse_discrete_quantile(
 384       dist,
 385       p,
 386       c,
 387       guess,
 388       multiplier,
 389       adder,
 390       tools::eps_tolerance<typename Dist::value_type>(policies::digits<typename Dist::value_type, typename Dist::policy_type>()),
 391       max_iter);
 392 }
 393
 394 template <class Dist>
 395 inline typename Dist::value_type
 396    inverse_discrete_quantile(
 397       const Dist& dist,
 398       const typename Dist::value_type& p,
 399       bool c,
 400       const typename Dist::value_type& guess,
 401       const typename Dist::value_type& multiplier,
 402       const typename Dist::value_type& adder,
 403       const policies::discrete_quantile<policies::integer_round_outwards>&,
 404       boost::uintmax_t& max_iter)
 405 {
 406    typedef typename Dist::value_type value_type;
 407    BOOST_MATH_STD_USING
 408    typename Dist::value_type pp = c ? 1 - p : p;
 409    if(pp <= pdf(dist, 0))
 410       return 0;
 411    //
 412    // What happens next depends on whether we're looking for an
 413    // upper or lower quantile:
 414    //
 415    if(pp < 0.5f)
 416       return round_to_floor(dist, do_inverse_discrete_quantile(
 417          dist,
 418          p,
 419          c,
 420          (guess < 1 ? value_type(1) : (value_type)floor(guess)),
 421          multiplier,
 422          adder,
 423          tools::equal_floor(),
 424          max_iter), p, c);
 425    // else:
 426    return round_to_ceil(dist, do_inverse_discrete_quantile(
 427       dist,
 428       p,
 429       c,
 430       (value_type)ceil(guess),
 431       multiplier,
 432       adder,
 433       tools::equal_ceil(),
 434       max_iter), p, c);
 435 }
 436
 437 template <class Dist>
 438 inline typename Dist::value_type
 439    inverse_discrete_quantile(
 440       const Dist& dist,
 441       const typename Dist::value_type& p,
 442       bool c,
 443       const typename Dist::value_type& guess,
 444       const typename Dist::value_type& multiplier,
 445       const typename Dist::value_type& adder,
 446       const policies::discrete_quantile<policies::integer_round_inwards>&,
 447       boost::uintmax_t& max_iter)
 448 {
 449    typedef typename Dist::value_type value_type;
 450    BOOST_MATH_STD_USING
 451    typename Dist::value_type pp = c ? 1 - p : p;
 452    if(pp <= pdf(dist, 0))
 453       return 0;
 454    //
 455    // What happens next depends on whether we're looking for an
 456    // upper or lower quantile:
 457    //
 458    if(pp < 0.5f)
 459       return round_to_ceil(dist, do_inverse_discrete_quantile(
 460          dist,
 461          p,
 462          c,
 463          ceil(guess),
 464          multiplier,
 465          adder,
 466          tools::equal_ceil(),
 467          max_iter), p, c);
 468    // else:
 469    return round_to_floor(dist, do_inverse_discrete_quantile(
 470       dist,
 471       p,
 472       c,
 473       (guess < 1 ? value_type(1) : floor(guess)),
 474       multiplier,
 475       adder,
 476       tools::equal_floor(),
 477       max_iter), p, c);
 478 }
 479
 480 template <class Dist>
 481 inline typename Dist::value_type
 482    inverse_discrete_quantile(
 483       const Dist& dist,
 484       const typename Dist::value_type& p,
 485       bool c,
 486       const typename Dist::value_type& guess,
 487       const typename Dist::value_type& multiplier,
 488       const typename Dist::value_type& adder,
 489       const policies::discrete_quantile<policies::integer_round_down>&,
 490       boost::uintmax_t& max_iter)
 491 {
 492    typedef typename Dist::value_type value_type;
 493    BOOST_MATH_STD_USING
 494    typename Dist::value_type pp = c ? 1 - p : p;
 495    if(pp <= pdf(dist, 0))
 496       return 0;
 497    return round_to_floor(dist, do_inverse_discrete_quantile(
 498       dist,
 499       p,
 500       c,
 501       (guess < 1 ? value_type(1) : floor(guess)),
 502       multiplier,
 503       adder,
 504       tools::equal_floor(),
 505       max_iter), p, c);
 506 }
 507
 508 template <class Dist>
 509 inline typename Dist::value_type
 510    inverse_discrete_quantile(
 511       const Dist& dist,
 512       const typename Dist::value_type& p,
 513       bool c,
 514       const typename Dist::value_type& guess,
 515       const typename Dist::value_type& multiplier,
 516       const typename Dist::value_type& adder,
 517       const policies::discrete_quantile<policies::integer_round_up>&,
 518       boost::uintmax_t& max_iter)
 519 {
 520    BOOST_MATH_STD_USING
 521    typename Dist::value_type pp = c ? 1 - p : p;
 522    if(pp <= pdf(dist, 0))
 523       return 0;
 524    return round_to_ceil(dist, do_inverse_discrete_quantile(
 525       dist,
 526       p,
 527       c,
 528       ceil(guess),
 529       multiplier,
 530       adder,
 531       tools::equal_ceil(),
 532       max_iter), p, c);
 533 }
 534
 535 template <class Dist>
 536 inline typename Dist::value_type
 537    inverse_discrete_quantile(
 538       const Dist& dist,
 539       const typename Dist::value_type& p,
 540       bool c,
 541       const typename Dist::value_type& guess,
 542       const typename Dist::value_type& multiplier,
 543       const typename Dist::value_type& adder,
 544       const policies::discrete_quantile<policies::integer_round_nearest>&,
 545       boost::uintmax_t& max_iter)
 546 {
 547    typedef typename Dist::value_type value_type;
 548    BOOST_MATH_STD_USING
 549    typename Dist::value_type pp = c ? 1 - p : p;
 550    if(pp <= pdf(dist, 0))
 551       return 0;
 552    //
 553    // Note that we adjust the guess to the nearest half-integer:
 554    // this increase the chances that we will bracket the root
 555    // with two results that both round to the same integer quickly.
 556    //
 557    return round_to_floor(dist, do_inverse_discrete_quantile(
 558       dist,
 559       p,
 560       c,
 561       (guess < 0.5f ? value_type(1.5f) : floor(guess + 0.5f) + 0.5f),
 562       multiplier,
 563       adder,
 564       tools::equal_nearest_integer(),
 565       max_iter) + 0.5f, p, c);
 566 }
 567
 568 }}} // namespaces
 569
 570 #endif // BOOST_MATH_DISTRIBUTIONS_DETAIL_INV_DISCRETE_QUANTILE
 571