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