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