diff options
Diffstat (limited to 'src/third_party/boost-1.69.0/boost/math/distributions/detail/inv_discrete_quantile.hpp')
| -rw-r--r-- | src/third_party/boost-1.69.0/boost/math/distributions/detail/inv_discrete_quantile.hpp | 571 |
1 files changed, 0 insertions, 571 deletions
diff --git a/src/third_party/boost-1.69.0/boost/math/distributions/detail/inv_discrete_quantile.hpp b/src/third_party/boost-1.69.0/boost/math/distributions/detail/inv_discrete_quantile.hpp deleted file mode 100644 index 23e00b8e036..00000000000 --- a/src/third_party/boost-1.69.0/boost/math/distributions/detail/inv_discrete_quantile.hpp +++ /dev/null @@ -1,571 +0,0 @@ -// Copyright John Maddock 2007. -// Use, modification and distribution are subject to the -// Boost Software License, Version 1.0. (See accompanying file -// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) - -#ifndef BOOST_MATH_DISTRIBUTIONS_DETAIL_INV_DISCRETE_QUANTILE -#define BOOST_MATH_DISTRIBUTIONS_DETAIL_INV_DISCRETE_QUANTILE - -#include <algorithm> - -namespace boost{ namespace math{ namespace detail{ - -// -// Functor for root finding algorithm: -// -template <class Dist> -struct distribution_quantile_finder -{ - typedef typename Dist::value_type value_type; - typedef typename Dist::policy_type policy_type; - - distribution_quantile_finder(const Dist d, value_type p, bool c) - : dist(d), target(p), comp(c) {} - - value_type operator()(value_type const& x) - { - return comp ? value_type(target - cdf(complement(dist, x))) : value_type(cdf(dist, x) - target); - } - -private: - Dist dist; - value_type target; - bool comp; -}; -// -// The purpose of adjust_bounds, is to toggle the last bit of the -// range so that both ends round to the same integer, if possible. -// If they do both round the same then we terminate the search -// for the root *very* quickly when finding an integer result. -// At the point that this function is called we know that "a" is -// below the root and "b" above it, so this change can not result -// in the root no longer being bracketed. -// -template <class Real, class Tol> -void adjust_bounds(Real& /* a */, Real& /* b */, Tol const& /* tol */){} - -template <class Real> -void adjust_bounds(Real& /* a */, Real& b, tools::equal_floor const& /* tol */) -{ - BOOST_MATH_STD_USING - b -= tools::epsilon<Real>() * b; -} - -template <class Real> -void adjust_bounds(Real& a, Real& /* b */, tools::equal_ceil const& /* tol */) -{ - BOOST_MATH_STD_USING - a += tools::epsilon<Real>() * a; -} - -template <class Real> -void adjust_bounds(Real& a, Real& b, tools::equal_nearest_integer const& /* tol */) -{ - BOOST_MATH_STD_USING - a += tools::epsilon<Real>() * a; - b -= tools::epsilon<Real>() * b; -} -// -// This is where all the work is done: -// -template <class Dist, class Tolerance> -typename Dist::value_type - do_inverse_discrete_quantile( - const Dist& dist, - const typename Dist::value_type& p, - bool comp, - typename Dist::value_type guess, - const typename Dist::value_type& multiplier, - typename Dist::value_type adder, - const Tolerance& tol, - boost::uintmax_t& max_iter) -{ - typedef typename Dist::value_type value_type; - typedef typename Dist::policy_type policy_type; - - static const char* function = "boost::math::do_inverse_discrete_quantile<%1%>"; - - BOOST_MATH_STD_USING - - distribution_quantile_finder<Dist> f(dist, p, comp); - // - // Max bounds of the distribution: - // - value_type min_bound, max_bound; - boost::math::tie(min_bound, max_bound) = support(dist); - - if(guess > max_bound) - guess = max_bound; - if(guess < min_bound) - guess = min_bound; - - value_type fa = f(guess); - boost::uintmax_t count = max_iter - 1; - value_type fb(fa), a(guess), b =0; // Compiler warning C4701: potentially uninitialized local variable 'b' used - - if(fa == 0) - return guess; - - // - // For small expected results, just use a linear search: - // - if(guess < 10) - { - b = a; - while((a < 10) && (fa * fb >= 0)) - { - if(fb <= 0) - { - a = b; - b = a + 1; - if(b > max_bound) - b = max_bound; - fb = f(b); - --count; - if(fb == 0) - return b; - if(a == b) - return b; // can't go any higher! - } - else - { - b = a; - a = (std::max)(value_type(b - 1), value_type(0)); - if(a < min_bound) - a = min_bound; - fa = f(a); - --count; - if(fa == 0) - return a; - if(a == b) - return a; // We can't go any lower than this! - } - } - } - // - // Try and bracket using a couple of additions first, - // we're assuming that "guess" is likely to be accurate - // to the nearest int or so: - // - else if(adder != 0) - { - // - // If we're looking for a large result, then bump "adder" up - // by a bit to increase our chances of bracketing the root: - // - //adder = (std::max)(adder, 0.001f * guess); - if(fa < 0) - { - b = a + adder; - if(b > max_bound) - b = max_bound; - } - else - { - b = (std::max)(value_type(a - adder), value_type(0)); - if(b < min_bound) - b = min_bound; - } - fb = f(b); - --count; - if(fb == 0) - return b; - if(count && (fa * fb >= 0)) - { - // - // We didn't bracket the root, try - // once more: - // - a = b; - fa = fb; - if(fa < 0) - { - b = a + adder; - if(b > max_bound) - b = max_bound; - } - else - { - b = (std::max)(value_type(a - adder), value_type(0)); - if(b < min_bound) - b = min_bound; - } - fb = f(b); - --count; - } - if(a > b) - { - using std::swap; - swap(a, b); - swap(fa, fb); - } - } - // - // If the root hasn't been bracketed yet, try again - // using the multiplier this time: - // - if((boost::math::sign)(fb) == (boost::math::sign)(fa)) - { - if(fa < 0) - { - // - // Zero is to the right of x2, so walk upwards - // until we find it: - // - while(((boost::math::sign)(fb) == (boost::math::sign)(fa)) && (a != b)) - { - if(count == 0) - return policies::raise_evaluation_error(function, "Unable to bracket root, last nearest value was %1%", b, policy_type()); - a = b; - fa = fb; - b *= multiplier; - if(b > max_bound) - b = max_bound; - fb = f(b); - --count; - BOOST_MATH_INSTRUMENT_CODE("a = " << a << " b = " << b << " fa = " << fa << " fb = " << fb << " count = " << count); - } - } - else - { - // - // Zero is to the left of a, so walk downwards - // until we find it: - // - while(((boost::math::sign)(fb) == (boost::math::sign)(fa)) && (a != b)) - { - if(fabs(a) < tools::min_value<value_type>()) - { - // Escape route just in case the answer is zero! - max_iter -= count; - max_iter += 1; - return 0; - } - if(count == 0) - return policies::raise_evaluation_error(function, "Unable to bracket root, last nearest value was %1%", a, policy_type()); - b = a; - fb = fa; - a /= multiplier; - if(a < min_bound) - a = min_bound; - fa = f(a); - --count; - BOOST_MATH_INSTRUMENT_CODE("a = " << a << " b = " << b << " fa = " << fa << " fb = " << fb << " count = " << count); - } - } - } - max_iter -= count; - if(fa == 0) - return a; - if(fb == 0) - return b; - if(a == b) - return b; // Ran out of bounds trying to bracket - there is no answer! - // - // Adjust bounds so that if we're looking for an integer - // result, then both ends round the same way: - // - adjust_bounds(a, b, tol); - // - // We don't want zero or denorm lower bounds: - // - if(a < tools::min_value<value_type>()) - a = tools::min_value<value_type>(); - // - // Go ahead and find the root: - // - std::pair<value_type, value_type> r = toms748_solve(f, a, b, fa, fb, tol, count, policy_type()); - max_iter += count; - BOOST_MATH_INSTRUMENT_CODE("max_iter = " << max_iter << " count = " << count); - return (r.first + r.second) / 2; -} -// -// Some special routine for rounding up and down: -// We want to check and see if we are very close to an integer, and if so test to see if -// that integer is an exact root of the cdf. We do this because our root finder only -// guarantees to find *a root*, and there can sometimes be many consecutive floating -// point values which are all roots. This is especially true if the target probability -// is very close 1. -// -template <class Dist> -inline typename Dist::value_type round_to_floor(const Dist& d, typename Dist::value_type result, typename Dist::value_type p, bool c) -{ - BOOST_MATH_STD_USING - typename Dist::value_type cc = ceil(result); - typename Dist::value_type pp = cc <= support(d).second ? c ? cdf(complement(d, cc)) : cdf(d, cc) : 1; - if(pp == p) - result = cc; - else - result = floor(result); - // - // Now find the smallest integer <= result for which we get an exact root: - // - while(result != 0) - { - cc = result - 1; - if(cc < support(d).first) - break; - pp = c ? cdf(complement(d, cc)) : cdf(d, cc); - if(pp == p) - result = cc; - else if(c ? pp > p : pp < p) - break; - result -= 1; - } - - return result; -} - -#ifdef BOOST_MSVC -#pragma warning(push) -#pragma warning(disable:4127) -#endif - -template <class Dist> -inline typename Dist::value_type round_to_ceil(const Dist& d, typename Dist::value_type result, typename Dist::value_type p, bool c) -{ - BOOST_MATH_STD_USING - typename Dist::value_type cc = floor(result); - typename Dist::value_type pp = cc >= support(d).first ? c ? cdf(complement(d, cc)) : cdf(d, cc) : 0; - if(pp == p) - result = cc; - else - result = ceil(result); - // - // Now find the largest integer >= result for which we get an exact root: - // - while(true) - { - cc = result + 1; - if(cc > support(d).second) - break; - pp = c ? cdf(complement(d, cc)) : cdf(d, cc); - if(pp == p) - result = cc; - else if(c ? pp < p : pp > p) - break; - result += 1; - } - - return result; -} - -#ifdef BOOST_MSVC -#pragma warning(pop) -#endif -// -// Now finally are the public API functions. -// There is one overload for each policy, -// each one is responsible for selecting the correct -// termination condition, and rounding the result -// to an int where required. -// -template <class Dist> -inline typename Dist::value_type - inverse_discrete_quantile( - const Dist& dist, - typename Dist::value_type p, - bool c, - const typename Dist::value_type& guess, - const typename Dist::value_type& multiplier, - const typename Dist::value_type& adder, - const policies::discrete_quantile<policies::real>&, - boost::uintmax_t& max_iter) -{ - if(p > 0.5) - { - p = 1 - p; - c = !c; - } - typename Dist::value_type pp = c ? 1 - p : p; - if(pp <= pdf(dist, 0)) - return 0; - return do_inverse_discrete_quantile( - dist, - p, - c, - guess, - multiplier, - adder, - tools::eps_tolerance<typename Dist::value_type>(policies::digits<typename Dist::value_type, typename Dist::policy_type>()), - max_iter); -} - -template <class Dist> -inline typename Dist::value_type - inverse_discrete_quantile( - const Dist& dist, - const typename Dist::value_type& p, - bool c, - const typename Dist::value_type& guess, - const typename Dist::value_type& multiplier, - const typename Dist::value_type& adder, - const policies::discrete_quantile<policies::integer_round_outwards>&, - boost::uintmax_t& max_iter) -{ - typedef typename Dist::value_type value_type; - BOOST_MATH_STD_USING - typename Dist::value_type pp = c ? 1 - p : p; - if(pp <= pdf(dist, 0)) - return 0; - // - // What happens next depends on whether we're looking for an - // upper or lower quantile: - // - if(pp < 0.5f) - return round_to_floor(dist, do_inverse_discrete_quantile( - dist, - p, - c, - (guess < 1 ? value_type(1) : (value_type)floor(guess)), - multiplier, - adder, - tools::equal_floor(), - max_iter), p, c); - // else: - return round_to_ceil(dist, do_inverse_discrete_quantile( - dist, - p, - c, - (value_type)ceil(guess), - multiplier, - adder, - tools::equal_ceil(), - max_iter), p, c); -} - -template <class Dist> -inline typename Dist::value_type - inverse_discrete_quantile( - const Dist& dist, - const typename Dist::value_type& p, - bool c, - const typename Dist::value_type& guess, - const typename Dist::value_type& multiplier, - const typename Dist::value_type& adder, - const policies::discrete_quantile<policies::integer_round_inwards>&, - boost::uintmax_t& max_iter) -{ - typedef typename Dist::value_type value_type; - BOOST_MATH_STD_USING - typename Dist::value_type pp = c ? 1 - p : p; - if(pp <= pdf(dist, 0)) - return 0; - // - // What happens next depends on whether we're looking for an - // upper or lower quantile: - // - if(pp < 0.5f) - return round_to_ceil(dist, do_inverse_discrete_quantile( - dist, - p, - c, - ceil(guess), - multiplier, - adder, - tools::equal_ceil(), - max_iter), p, c); - // else: - return round_to_floor(dist, do_inverse_discrete_quantile( - dist, - p, - c, - (guess < 1 ? value_type(1) : floor(guess)), - multiplier, - adder, - tools::equal_floor(), - max_iter), p, c); -} - -template <class Dist> -inline typename Dist::value_type - inverse_discrete_quantile( - const Dist& dist, - const typename Dist::value_type& p, - bool c, - const typename Dist::value_type& guess, - const typename Dist::value_type& multiplier, - const typename Dist::value_type& adder, - const policies::discrete_quantile<policies::integer_round_down>&, - boost::uintmax_t& max_iter) -{ - typedef typename Dist::value_type value_type; - BOOST_MATH_STD_USING - typename Dist::value_type pp = c ? 1 - p : p; - if(pp <= pdf(dist, 0)) - return 0; - return round_to_floor(dist, do_inverse_discrete_quantile( - dist, - p, - c, - (guess < 1 ? value_type(1) : floor(guess)), - multiplier, - adder, - tools::equal_floor(), - max_iter), p, c); -} - -template <class Dist> -inline typename Dist::value_type - inverse_discrete_quantile( - const Dist& dist, - const typename Dist::value_type& p, - bool c, - const typename Dist::value_type& guess, - const typename Dist::value_type& multiplier, - const typename Dist::value_type& adder, - const policies::discrete_quantile<policies::integer_round_up>&, - boost::uintmax_t& max_iter) -{ - BOOST_MATH_STD_USING - typename Dist::value_type pp = c ? 1 - p : p; - if(pp <= pdf(dist, 0)) - return 0; - return round_to_ceil(dist, do_inverse_discrete_quantile( - dist, - p, - c, - ceil(guess), - multiplier, - adder, - tools::equal_ceil(), - max_iter), p, c); -} - -template <class Dist> -inline typename Dist::value_type - inverse_discrete_quantile( - const Dist& dist, - const typename Dist::value_type& p, - bool c, - const typename Dist::value_type& guess, - const typename Dist::value_type& multiplier, - const typename Dist::value_type& adder, - const policies::discrete_quantile<policies::integer_round_nearest>&, - boost::uintmax_t& max_iter) -{ - typedef typename Dist::value_type value_type; - BOOST_MATH_STD_USING - typename Dist::value_type pp = c ? 1 - p : p; - if(pp <= pdf(dist, 0)) - return 0; - // - // Note that we adjust the guess to the nearest half-integer: - // this increase the chances that we will bracket the root - // with two results that both round to the same integer quickly. - // - return round_to_floor(dist, do_inverse_discrete_quantile( - dist, - p, - c, - (guess < 0.5f ? value_type(1.5f) : floor(guess + 0.5f) + 0.5f), - multiplier, - adder, - tools::equal_nearest_integer(), - max_iter) + 0.5f, p, c); -} - -}}} // namespaces - -#endif // BOOST_MATH_DISTRIBUTIONS_DETAIL_INV_DISCRETE_QUANTILE - |