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