Boost logo

Boost-Commit :

Subject: [Boost-commit] svn:boost r86205 - in trunk/boost/math/distributions: . detail
From: john_at_[hidden]
Date: 2013-10-08 13:17:28


Author: johnmaddock
Date: 2013-10-08 13:17:27 EDT (Tue, 08 Oct 2013)
New Revision: 86205
URL: http://svn.boost.org/trac/boost/changeset/86205

Log:
Improve quantiles of discrete distributions to round trip integers more often.
Fixes #9183.

Text files modified:
   trunk/boost/math/distributions/binomial.hpp | 10 +-
   trunk/boost/math/distributions/detail/inv_discrete_quantile.hpp | 158 +++++++++++++++++++++++++++++----------
   trunk/boost/math/distributions/negative_binomial.hpp | 4
   trunk/boost/math/distributions/poisson.hpp | 66 ----------------
   4 files changed, 125 insertions(+), 113 deletions(-)

Modified: trunk/boost/math/distributions/binomial.hpp
==============================================================================
--- trunk/boost/math/distributions/binomial.hpp Tue Oct 8 13:15:17 2013 (r86204)
+++ trunk/boost/math/distributions/binomial.hpp 2013-10-08 13:17:27 EDT (Tue, 08 Oct 2013) (r86205)
@@ -196,7 +196,7 @@
          }
 
       template <class RealType, class Policy>
- RealType quantile_imp(const binomial_distribution<RealType, Policy>& dist, const RealType& p, const RealType& q)
+ RealType quantile_imp(const binomial_distribution<RealType, Policy>& dist, const RealType& p, const RealType& q, bool comp)
       { // Quantile or Percent Point Binomial function.
         // Return the number of expected successes k,
         // for a given probability p.
@@ -264,8 +264,8 @@
         boost::uintmax_t max_iter = policies::get_max_root_iterations<Policy>();
         return detail::inverse_discrete_quantile(
             dist,
- p,
- q,
+ comp ? q : p,
+ comp,
             guess,
             factor,
             RealType(1),
@@ -653,13 +653,13 @@
       template <class RealType, class Policy>
       inline RealType quantile(const binomial_distribution<RealType, Policy>& dist, const RealType& p)
       {
- return binomial_detail::quantile_imp(dist, p, RealType(1-p));
+ return binomial_detail::quantile_imp(dist, p, RealType(1-p), false);
       } // quantile
 
       template <class RealType, class Policy>
       RealType quantile(const complemented2_type<binomial_distribution<RealType, Policy>, RealType>& c)
       {
- return binomial_detail::quantile_imp(c.dist, RealType(1-c.param), c.param);
+ return binomial_detail::quantile_imp(c.dist, RealType(1-c.param), c.param, true);
       } // quantile
 
       template <class RealType, class Policy>

Modified: trunk/boost/math/distributions/detail/inv_discrete_quantile.hpp
==============================================================================
--- trunk/boost/math/distributions/detail/inv_discrete_quantile.hpp Tue Oct 8 13:15:17 2013 (r86204)
+++ trunk/boost/math/distributions/detail/inv_discrete_quantile.hpp 2013-10-08 13:17:27 EDT (Tue, 08 Oct 2013) (r86205)
@@ -19,8 +19,8 @@
    typedef typename Dist::value_type value_type;
    typedef typename Dist::policy_type policy_type;
 
- distribution_quantile_finder(const Dist d, value_type p, value_type q)
- : dist(d), target(p < q ? p : q), comp(p < q ? false : true) {}
+ distribution_quantile_finder(const Dist d, value_type p, bool c)
+ : dist(d), target(p), comp(c) {}
 
    value_type operator()(value_type const& x)
    {
@@ -73,7 +73,7 @@
    do_inverse_discrete_quantile(
       const Dist& dist,
       const typename Dist::value_type& p,
- const typename Dist::value_type& q,
+ bool comp,
       typename Dist::value_type guess,
       const typename Dist::value_type& multiplier,
       typename Dist::value_type adder,
@@ -87,7 +87,7 @@
 
    BOOST_MATH_STD_USING
 
- distribution_quantile_finder<Dist> f(dist, p, q);
+ distribution_quantile_finder<Dist> f(dist, p, comp);
    //
    // Max bounds of the distribution:
    //
@@ -280,6 +280,70 @@
    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;
+ typename Dist::value_type 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;
+}
+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;
+ typename Dist::value_type 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;
+}
+//
 // Now finally are the public API functions.
 // There is one overload for each policy,
 // each one is responsible for selecting the correct
@@ -290,20 +354,26 @@
 inline typename Dist::value_type
    inverse_discrete_quantile(
       const Dist& dist,
- const typename Dist::value_type& p,
- const typename Dist::value_type& q,
+ 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 <= pdf(dist, 0))
+ 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,
- q,
+ c,
       guess,
       multiplier,
       adder,
@@ -316,7 +386,7 @@
    inverse_discrete_quantile(
       const Dist& dist,
       const typename Dist::value_type& p,
- const typename Dist::value_type& q,
+ bool c,
       const typename Dist::value_type& guess,
       const typename Dist::value_type& multiplier,
       const typename Dist::value_type& adder,
@@ -325,32 +395,33 @@
 {
    typedef typename Dist::value_type value_type;
    BOOST_MATH_STD_USING
- if(p <= pdf(dist, 0))
+ 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(p < 0.5f)
- return floor(do_inverse_discrete_quantile(
+ if(pp < 0.5f)
+ return round_to_floor(dist, do_inverse_discrete_quantile(
          dist,
          p,
- q,
+ c,
          (guess < 1 ? value_type(1) : (value_type)floor(guess)),
          multiplier,
          adder,
          tools::equal_floor(),
- max_iter));
+ max_iter), p, c);
    // else:
- return ceil(do_inverse_discrete_quantile(
+ return round_to_ceil(dist, do_inverse_discrete_quantile(
       dist,
       p,
- q,
+ c,
       (value_type)ceil(guess),
       multiplier,
       adder,
       tools::equal_ceil(),
- max_iter));
+ max_iter), p, c);
 }
 
 template <class Dist>
@@ -358,7 +429,7 @@
    inverse_discrete_quantile(
       const Dist& dist,
       const typename Dist::value_type& p,
- const typename Dist::value_type& q,
+ bool c,
       const typename Dist::value_type& guess,
       const typename Dist::value_type& multiplier,
       const typename Dist::value_type& adder,
@@ -367,32 +438,33 @@
 {
    typedef typename Dist::value_type value_type;
    BOOST_MATH_STD_USING
- if(p <= pdf(dist, 0))
+ 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(p < 0.5f)
- return ceil(do_inverse_discrete_quantile(
+ if(pp < 0.5f)
+ return round_to_ceil(dist, do_inverse_discrete_quantile(
          dist,
          p,
- q,
+ c,
          ceil(guess),
          multiplier,
          adder,
          tools::equal_ceil(),
- max_iter));
+ max_iter), p, c);
    // else:
- return floor(do_inverse_discrete_quantile(
+ return round_to_floor(dist, do_inverse_discrete_quantile(
       dist,
       p,
- q,
+ c,
       (guess < 1 ? value_type(1) : floor(guess)),
       multiplier,
       adder,
       tools::equal_floor(),
- max_iter));
+ max_iter), p, c);
 }
 
 template <class Dist>
@@ -400,7 +472,7 @@
    inverse_discrete_quantile(
       const Dist& dist,
       const typename Dist::value_type& p,
- const typename Dist::value_type& q,
+ bool c,
       const typename Dist::value_type& guess,
       const typename Dist::value_type& multiplier,
       const typename Dist::value_type& adder,
@@ -409,17 +481,18 @@
 {
    typedef typename Dist::value_type value_type;
    BOOST_MATH_STD_USING
- if(p <= pdf(dist, 0))
+ typename Dist::value_type pp = c ? 1 - p : p;
+ if(pp <= pdf(dist, 0))
       return 0;
- return floor(do_inverse_discrete_quantile(
+ return round_to_floor(dist, do_inverse_discrete_quantile(
       dist,
       p,
- q,
+ c,
       (guess < 1 ? value_type(1) : floor(guess)),
       multiplier,
       adder,
       tools::equal_floor(),
- max_iter));
+ max_iter), p, c);
 }
 
 template <class Dist>
@@ -427,7 +500,7 @@
    inverse_discrete_quantile(
       const Dist& dist,
       const typename Dist::value_type& p,
- const typename Dist::value_type& q,
+ bool c,
       const typename Dist::value_type& guess,
       const typename Dist::value_type& multiplier,
       const typename Dist::value_type& adder,
@@ -435,17 +508,18 @@
       boost::uintmax_t& max_iter)
 {
    BOOST_MATH_STD_USING
- if(p <= pdf(dist, 0))
+ typename Dist::value_type pp = c ? 1 - p : p;
+ if(pp <= pdf(dist, 0))
       return 0;
- return ceil(do_inverse_discrete_quantile(
+ return round_to_ceil(dist, do_inverse_discrete_quantile(
       dist,
       p,
- q,
+ c,
       ceil(guess),
       multiplier,
       adder,
       tools::equal_ceil(),
- max_iter));
+ max_iter), p, c);
 }
 
 template <class Dist>
@@ -453,7 +527,7 @@
    inverse_discrete_quantile(
       const Dist& dist,
       const typename Dist::value_type& p,
- const typename Dist::value_type& q,
+ bool c,
       const typename Dist::value_type& guess,
       const typename Dist::value_type& multiplier,
       const typename Dist::value_type& adder,
@@ -462,26 +536,26 @@
 {
    typedef typename Dist::value_type value_type;
    BOOST_MATH_STD_USING
- if(p <= pdf(dist, 0))
+ 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 floor(do_inverse_discrete_quantile(
+ return round_to_floor(dist, do_inverse_discrete_quantile(
       dist,
       p,
- q,
+ c,
       (guess < 0.5f ? value_type(1.5f) : floor(guess + 0.5f) + 0.5f),
       multiplier,
       adder,
       tools::equal_nearest_integer(),
- max_iter) + 0.5f);
+ max_iter) + 0.5f, p, c);
 }
 
 }}} // namespaces
 
 #endif // BOOST_MATH_DISTRIBUTIONS_DETAIL_INV_DISCRETE_QUANTILE
 
-

Modified: trunk/boost/math/distributions/negative_binomial.hpp
==============================================================================
--- trunk/boost/math/distributions/negative_binomial.hpp Tue Oct 8 13:15:17 2013 (r86204)
+++ trunk/boost/math/distributions/negative_binomial.hpp 2013-10-08 13:17:27 EDT (Tue, 08 Oct 2013) (r86205)
@@ -488,7 +488,7 @@
       return detail::inverse_discrete_quantile(
          dist,
          P,
- 1-P,
+ false,
          guess,
          factor,
          RealType(1),
@@ -564,8 +564,8 @@
        typedef typename Policy::discrete_quantile_type discrete_type;
        return detail::inverse_discrete_quantile(
           dist,
- 1-Q,
           Q,
+ true,
           guess,
           factor,
           RealType(1),

Modified: trunk/boost/math/distributions/poisson.hpp
==============================================================================
--- trunk/boost/math/distributions/poisson.hpp Tue Oct 8 13:15:17 2013 (r86204)
+++ trunk/boost/math/distributions/poisson.hpp 2013-10-08 13:17:27 EDT (Tue, 08 Oct 2013) (r86205)
@@ -52,68 +52,6 @@
 {
   namespace math
   {
- namespace detail{
- template <class Dist>
- inline typename Dist::value_type
- inverse_discrete_quantile(
- const Dist& dist,
- const typename Dist::value_type& p,
- 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);
- template <class Dist>
- inline typename Dist::value_type
- inverse_discrete_quantile(
- const Dist& dist,
- const typename Dist::value_type& p,
- 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);
- template <class Dist>
- inline typename Dist::value_type
- inverse_discrete_quantile(
- const Dist& dist,
- const typename Dist::value_type& p,
- 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);
- template <class Dist>
- inline typename Dist::value_type
- inverse_discrete_quantile(
- const Dist& dist,
- const typename Dist::value_type& p,
- 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);
- template <class Dist>
- inline typename Dist::value_type
- inverse_discrete_quantile(
- const Dist& dist,
- const typename Dist::value_type& p,
- 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);
- template <class Dist>
- inline typename Dist::value_type
- inverse_discrete_quantile(
- const Dist& dist,
- const typename Dist::value_type& p,
- 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);
- }
     namespace poisson_detail
     {
       // Common error checking routines for Poisson distribution functions.
@@ -496,7 +434,7 @@
       return detail::inverse_discrete_quantile(
          dist,
          p,
- 1-p,
+ false,
          guess,
          factor,
          RealType(1),
@@ -565,8 +503,8 @@
 
       return detail::inverse_discrete_quantile(
          dist,
- 1-q,
          q,
+ true,
          guess,
          factor,
          RealType(1),


Boost-Commit list run by bdawes at acm.org, david.abrahams at rcn.com, gregod at cs.rpi.edu, cpdaniel at pacbell.net, john at johnmaddock.co.uk