| // 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 |
| |