Imported existing code
This commit is contained in:
Hazim Gazov
@@ -0,0 +1,157 @@
|
||||
// Copyright John Maddock 2006, 2007.
|
||||
// Copyright Paul A. Bristow 2006, 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_COMMON_ERROR_HANDLING_HPP
|
||||
#define BOOST_MATH_DISTRIBUTIONS_COMMON_ERROR_HANDLING_HPP
|
||||
|
||||
#include <boost/math/policies/error_handling.hpp>
|
||||
#include <boost/math/special_functions/fpclassify.hpp>
|
||||
// using boost::math::isfinite;
|
||||
|
||||
namespace boost{ namespace math{ namespace detail
|
||||
{
|
||||
|
||||
template <class RealType, class Policy>
|
||||
inline bool check_probability(const char* function, RealType const& prob, RealType* result, const Policy& pol)
|
||||
{
|
||||
if((prob < 0) || (prob > 1) || !(boost::math::isfinite)(prob))
|
||||
{
|
||||
*result = policies::raise_domain_error<RealType>(
|
||||
function,
|
||||
"Probability argument is %1%, but must be >= 0 and <= 1 !", prob, pol);
|
||||
return false;
|
||||
}
|
||||
return true;
|
||||
}
|
||||
|
||||
template <class RealType, class Policy>
|
||||
inline bool check_df(const char* function, RealType const& df, RealType* result, const Policy& pol)
|
||||
{
|
||||
if((df <= 0) || !(boost::math::isfinite)(df))
|
||||
{
|
||||
*result = policies::raise_domain_error<RealType>(
|
||||
function,
|
||||
"Degrees of freedom argument is %1%, but must be > 0 !", df, pol);
|
||||
return false;
|
||||
}
|
||||
return true;
|
||||
}
|
||||
|
||||
template <class RealType, class Policy>
|
||||
inline bool check_scale(
|
||||
const char* function,
|
||||
RealType scale,
|
||||
RealType* result,
|
||||
const Policy& pol)
|
||||
{
|
||||
if((scale <= 0) || !(boost::math::isfinite)(scale))
|
||||
{ // Assume scale == 0 is NOT valid for any distribution.
|
||||
*result = policies::raise_domain_error<RealType>(
|
||||
function,
|
||||
"Scale parameter is %1%, but must be > 0 !", scale, pol);
|
||||
return false;
|
||||
}
|
||||
return true;
|
||||
}
|
||||
|
||||
template <class RealType, class Policy>
|
||||
inline bool check_location(
|
||||
const char* function,
|
||||
RealType location,
|
||||
RealType* result,
|
||||
const Policy& pol)
|
||||
{
|
||||
if(!(boost::math::isfinite)(location))
|
||||
{
|
||||
*result = policies::raise_domain_error<RealType>(
|
||||
function,
|
||||
"Location parameter is %1%, but must be finite!", location, pol);
|
||||
return false;
|
||||
}
|
||||
return true;
|
||||
}
|
||||
|
||||
template <class RealType, class Policy>
|
||||
inline bool check_x(
|
||||
const char* function,
|
||||
RealType x,
|
||||
RealType* result,
|
||||
const Policy& pol)
|
||||
{
|
||||
if(!(boost::math::isfinite)(x))
|
||||
{
|
||||
*result = policies::raise_domain_error<RealType>(
|
||||
function,
|
||||
"Random variate x is %1%, but must be finite!", x, pol);
|
||||
return false;
|
||||
}
|
||||
return true;
|
||||
// Note that this test catches both infinity and NaN.
|
||||
// Some special cases permit x to be infinite, so these must be tested 1st,
|
||||
// leaving this test to catch any NaNs. see Normal and cauchy for example.
|
||||
}
|
||||
|
||||
template <class RealType, class Policy>
|
||||
inline bool check_positive_x(
|
||||
const char* function,
|
||||
RealType x,
|
||||
RealType* result,
|
||||
const Policy& pol)
|
||||
{
|
||||
if(!(boost::math::isfinite)(x) || (x < 0))
|
||||
{
|
||||
*result = policies::raise_domain_error<RealType>(
|
||||
function,
|
||||
"Random variate x is %1%, but must be finite and >= 0!", x, pol);
|
||||
return false;
|
||||
}
|
||||
return true;
|
||||
// Note that this test catches both infinity and NaN.
|
||||
// Some special cases permit x to be infinite, so these must be tested 1st,
|
||||
// leaving this test to catch any NaNs. see Normal and cauchy for example.
|
||||
}
|
||||
|
||||
template <class RealType, class Policy>
|
||||
inline bool check_non_centrality(
|
||||
const char* function,
|
||||
RealType ncp,
|
||||
RealType* result,
|
||||
const Policy& pol)
|
||||
{
|
||||
if((ncp < 0) || !(boost::math::isfinite)(ncp))
|
||||
{ // Assume scale == 0 is NOT valid for any distribution.
|
||||
*result = policies::raise_domain_error<RealType>(
|
||||
function,
|
||||
"Non centrality parameter is %1%, but must be > 0 !", ncp, pol);
|
||||
return false;
|
||||
}
|
||||
return true;
|
||||
}
|
||||
|
||||
template <class RealType, class Policy>
|
||||
inline bool check_finite(
|
||||
const char* function,
|
||||
RealType x,
|
||||
RealType* result,
|
||||
const Policy& pol)
|
||||
{
|
||||
if(!(boost::math::isfinite)(x))
|
||||
{ // Assume scale == 0 is NOT valid for any distribution.
|
||||
*result = policies::raise_domain_error<RealType>(
|
||||
function,
|
||||
"Parameter is %1%, but must be finite !", x, pol);
|
||||
return false;
|
||||
}
|
||||
return true;
|
||||
}
|
||||
|
||||
} // namespace detail
|
||||
} // namespace math
|
||||
} // namespace boost
|
||||
|
||||
#endif // BOOST_MATH_DISTRIBUTIONS_COMMON_ERROR_HANDLING_HPP
|
||||
@@ -0,0 +1,163 @@
|
||||
// Copyright John Maddock 2006.
|
||||
// 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_STATS_DERIVED_HPP
|
||||
#define BOOST_STATS_DERIVED_HPP
|
||||
|
||||
// This file implements various common properties of distributions
|
||||
// that can be implemented in terms of other properties:
|
||||
// variance OR standard deviation (see note below),
|
||||
// hazard, cumulative hazard (chf), coefficient_of_variation.
|
||||
//
|
||||
// Note that while both variance and standard_deviation are provided
|
||||
// here, each distribution MUST SPECIALIZE AT LEAST ONE OF THESE
|
||||
// otherwise these two versions will just call each other over and over
|
||||
// until stack space runs out ...
|
||||
|
||||
// Of course there may be more efficient means of implementing these
|
||||
// that are specific to a particular distribution, but these generic
|
||||
// versions give these properties "for free" with most distributions.
|
||||
//
|
||||
// In order to make use of this header, it must be included AT THE END
|
||||
// of the distribution header, AFTER the distribution and its core
|
||||
// property accessors have been defined: this is so that compilers
|
||||
// that implement 2-phase lookup and early-type-checking of templates
|
||||
// can find the definitions refered to herein.
|
||||
//
|
||||
|
||||
#include <boost/type_traits/is_same.hpp>
|
||||
#include <boost/static_assert.hpp>
|
||||
|
||||
#ifdef BOOST_MSVC
|
||||
# pragma warning(push)
|
||||
# pragma warning(disable: 4723) // potential divide by 0
|
||||
// Suppressing spurious warning in coefficient_of_variation
|
||||
#endif
|
||||
|
||||
namespace boost{ namespace math{
|
||||
|
||||
template <class Distribution>
|
||||
typename Distribution::value_type variance(const Distribution& dist);
|
||||
|
||||
template <class Distribution>
|
||||
inline typename Distribution::value_type standard_deviation(const Distribution& dist)
|
||||
{
|
||||
BOOST_MATH_STD_USING // ADL of sqrt.
|
||||
return sqrt(variance(dist));
|
||||
}
|
||||
|
||||
template <class Distribution>
|
||||
inline typename Distribution::value_type variance(const Distribution& dist)
|
||||
{
|
||||
typename Distribution::value_type result = standard_deviation(dist);
|
||||
return result * result;
|
||||
}
|
||||
|
||||
template <class Distribution, class RealType>
|
||||
inline typename Distribution::value_type hazard(const Distribution& dist, const RealType& x)
|
||||
{ // hazard function
|
||||
// http://www.itl.nist.gov/div898/handbook/eda/section3/eda362.htm#HAZ
|
||||
typedef typename Distribution::value_type value_type;
|
||||
typedef typename Distribution::policy_type policy_type;
|
||||
value_type p = cdf(complement(dist, x));
|
||||
value_type d = pdf(dist, x);
|
||||
if(d > p * tools::max_value<value_type>())
|
||||
return policies::raise_overflow_error<value_type>(
|
||||
"boost::math::hazard(const Distribution&, %1%)", 0, policy_type());
|
||||
if(d == 0)
|
||||
{
|
||||
// This protects against 0/0, but is it the right thing to do?
|
||||
return 0;
|
||||
}
|
||||
return d / p;
|
||||
}
|
||||
|
||||
template <class Distribution, class RealType>
|
||||
inline typename Distribution::value_type chf(const Distribution& dist, const RealType& x)
|
||||
{ // cumulative hazard function.
|
||||
// http://www.itl.nist.gov/div898/handbook/eda/section3/eda362.htm#HAZ
|
||||
BOOST_MATH_STD_USING
|
||||
return -log(cdf(complement(dist, x)));
|
||||
}
|
||||
|
||||
template <class Distribution>
|
||||
inline typename Distribution::value_type coefficient_of_variation(const Distribution& dist)
|
||||
{
|
||||
typedef typename Distribution::value_type value_type;
|
||||
typedef typename Distribution::policy_type policy_type;
|
||||
|
||||
using std::abs;
|
||||
|
||||
value_type m = mean(dist);
|
||||
value_type d = standard_deviation(dist);
|
||||
if((abs(m) < 1) && (d > abs(m) * tools::max_value<value_type>()))
|
||||
{ // Checks too that m is not zero,
|
||||
return policies::raise_overflow_error<value_type>("boost::math::coefficient_of_variation(const Distribution&, %1%)", 0, policy_type());
|
||||
}
|
||||
return d / m; // so MSVC warning on zerodivide is spurious, and suppressed.
|
||||
}
|
||||
//
|
||||
// Next follow overloads of some of the standard accessors with mixed
|
||||
// argument types. We just use a typecast to forward on to the "real"
|
||||
// implementation with all arguments of the same type:
|
||||
//
|
||||
template <class Distribution, class RealType>
|
||||
inline typename Distribution::value_type pdf(const Distribution& dist, const RealType& x)
|
||||
{
|
||||
typedef typename Distribution::value_type value_type;
|
||||
return pdf(dist, static_cast<value_type>(x));
|
||||
}
|
||||
template <class Distribution, class RealType>
|
||||
inline typename Distribution::value_type cdf(const Distribution& dist, const RealType& x)
|
||||
{
|
||||
typedef typename Distribution::value_type value_type;
|
||||
return cdf(dist, static_cast<value_type>(x));
|
||||
}
|
||||
template <class Distribution, class RealType>
|
||||
inline typename Distribution::value_type quantile(const Distribution& dist, const RealType& x)
|
||||
{
|
||||
typedef typename Distribution::value_type value_type;
|
||||
return quantile(dist, static_cast<value_type>(x));
|
||||
}
|
||||
/*
|
||||
template <class Distribution, class RealType>
|
||||
inline typename Distribution::value_type chf(const Distribution& dist, const RealType& x)
|
||||
{
|
||||
typedef typename Distribution::value_type value_type;
|
||||
return chf(dist, static_cast<value_type>(x));
|
||||
}
|
||||
*/
|
||||
template <class Distribution, class RealType>
|
||||
inline typename Distribution::value_type cdf(const complemented2_type<Distribution, RealType>& c)
|
||||
{
|
||||
typedef typename Distribution::value_type value_type;
|
||||
return cdf(complement(c.dist, static_cast<value_type>(c.param)));
|
||||
}
|
||||
|
||||
template <class Distribution, class RealType>
|
||||
inline typename Distribution::value_type quantile(const complemented2_type<Distribution, RealType>& c)
|
||||
{
|
||||
typedef typename Distribution::value_type value_type;
|
||||
return quantile(complement(c.dist, static_cast<value_type>(c.param)));
|
||||
}
|
||||
|
||||
template <class Dist>
|
||||
inline typename Dist::value_type median(const Dist& d)
|
||||
{ // median - default definition for those distributions for which a
|
||||
// simple closed form is not known,
|
||||
// and for which a domain_error and/or NaN generating function is NOT defined.
|
||||
typedef typename Dist::value_type value_type;
|
||||
return quantile(d, static_cast<value_type>(0.5f));
|
||||
}
|
||||
|
||||
} // namespace math
|
||||
} // namespace boost
|
||||
|
||||
|
||||
#ifdef BOOST_MSVC
|
||||
# pragma warning(pop)
|
||||
#endif
|
||||
|
||||
#endif // BOOST_STATS_DERIVED_HPP
|
||||
@@ -0,0 +1,149 @@
|
||||
// Copyright John Maddock 2008.
|
||||
|
||||
// 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_MODE_HPP
|
||||
#define BOOST_MATH_DISTRIBUTIONS_DETAIL_MODE_HPP
|
||||
|
||||
#include <boost/math/tools/minima.hpp> // function minimization for mode
|
||||
#include <boost/math/policies/error_handling.hpp>
|
||||
#include <boost/math/distributions/fwd.hpp>
|
||||
|
||||
namespace boost{ namespace math{ namespace detail{
|
||||
|
||||
template <class Dist>
|
||||
struct pdf_minimizer
|
||||
{
|
||||
pdf_minimizer(const Dist& d)
|
||||
: dist(d) {}
|
||||
|
||||
typename Dist::value_type operator()(const typename Dist::value_type& x)
|
||||
{
|
||||
return -pdf(dist, x);
|
||||
}
|
||||
private:
|
||||
Dist dist;
|
||||
};
|
||||
|
||||
template <class Dist>
|
||||
typename Dist::value_type generic_find_mode(const Dist& dist, typename Dist::value_type guess, const char* function, typename Dist::value_type step = 0)
|
||||
{
|
||||
BOOST_MATH_STD_USING
|
||||
typedef typename Dist::value_type value_type;
|
||||
typedef typename Dist::policy_type policy_type;
|
||||
//
|
||||
// Need to begin by bracketing the maxima of the PDF:
|
||||
//
|
||||
value_type maxval;
|
||||
value_type upper_bound = guess;
|
||||
value_type lower_bound;
|
||||
value_type v = pdf(dist, guess);
|
||||
if(v == 0)
|
||||
{
|
||||
//
|
||||
// Oops we don't know how to handle this, or even in which
|
||||
// direction we should move in, treat as an evaluation error:
|
||||
//
|
||||
policies::raise_evaluation_error(
|
||||
function,
|
||||
"Could not locate a starting location for the search for the mode, original guess was %1%", guess, policy_type());
|
||||
}
|
||||
do
|
||||
{
|
||||
maxval = v;
|
||||
if(step != 0)
|
||||
upper_bound += step;
|
||||
else
|
||||
upper_bound *= 2;
|
||||
v = pdf(dist, upper_bound);
|
||||
}while(maxval < v);
|
||||
|
||||
lower_bound = upper_bound;
|
||||
do
|
||||
{
|
||||
maxval = v;
|
||||
if(step != 0)
|
||||
lower_bound -= step;
|
||||
else
|
||||
lower_bound /= 2;
|
||||
v = pdf(dist, lower_bound);
|
||||
}while(maxval < v);
|
||||
|
||||
boost::uintmax_t max_iter = policies::get_max_root_iterations<policy_type>();
|
||||
|
||||
value_type result = tools::brent_find_minima(
|
||||
pdf_minimizer<Dist>(dist),
|
||||
lower_bound,
|
||||
upper_bound,
|
||||
policies::digits<value_type, policy_type>(),
|
||||
max_iter).first;
|
||||
if(max_iter >= policies::get_max_root_iterations<policy_type>())
|
||||
{
|
||||
return policies::raise_evaluation_error<value_type>(
|
||||
function,
|
||||
"Unable to locate solution in a reasonable time:"
|
||||
" either there is no answer to the mode of the distribution"
|
||||
" or the answer is infinite. Current best guess is %1%", result, policy_type());
|
||||
}
|
||||
return result;
|
||||
}
|
||||
//
|
||||
// As above,but confined to the interval [0,1]:
|
||||
//
|
||||
template <class Dist>
|
||||
typename Dist::value_type generic_find_mode_01(const Dist& dist, typename Dist::value_type guess, const char* function)
|
||||
{
|
||||
BOOST_MATH_STD_USING
|
||||
typedef typename Dist::value_type value_type;
|
||||
typedef typename Dist::policy_type policy_type;
|
||||
//
|
||||
// Need to begin by bracketing the maxima of the PDF:
|
||||
//
|
||||
value_type maxval;
|
||||
value_type upper_bound = guess;
|
||||
value_type lower_bound;
|
||||
value_type v = pdf(dist, guess);
|
||||
do
|
||||
{
|
||||
maxval = v;
|
||||
upper_bound = 1 - (1 - upper_bound) / 2;
|
||||
if(upper_bound == 1)
|
||||
return 1;
|
||||
v = pdf(dist, upper_bound);
|
||||
}while(maxval < v);
|
||||
|
||||
lower_bound = upper_bound;
|
||||
do
|
||||
{
|
||||
maxval = v;
|
||||
lower_bound /= 2;
|
||||
if(lower_bound < tools::min_value<value_type>())
|
||||
return 0;
|
||||
v = pdf(dist, lower_bound);
|
||||
}while(maxval < v);
|
||||
|
||||
boost::uintmax_t max_iter = policies::get_max_root_iterations<policy_type>();
|
||||
|
||||
value_type result = tools::brent_find_minima(
|
||||
pdf_minimizer<Dist>(dist),
|
||||
lower_bound,
|
||||
upper_bound,
|
||||
policies::digits<value_type, policy_type>(),
|
||||
max_iter).first;
|
||||
if(max_iter >= policies::get_max_root_iterations<policy_type>())
|
||||
{
|
||||
return policies::raise_evaluation_error<value_type>(
|
||||
function,
|
||||
"Unable to locate solution in a reasonable time:"
|
||||
" either there is no answer to the mode of the distribution"
|
||||
" or the answer is infinite. Current best guess is %1%", result, policy_type());
|
||||
}
|
||||
return result;
|
||||
}
|
||||
|
||||
}}} // namespaces
|
||||
|
||||
#endif // BOOST_MATH_DISTRIBUTIONS_DETAIL_MODE_HPP
|
||||
@@ -0,0 +1,91 @@
|
||||
// Copyright John Maddock 2008.
|
||||
// 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_DISTIBUTIONS_DETAIL_GENERIC_QUANTILE_HPP
|
||||
#define BOOST_MATH_DISTIBUTIONS_DETAIL_GENERIC_QUANTILE_HPP
|
||||
|
||||
namespace boost{ namespace math{ namespace detail{
|
||||
|
||||
template <class Dist>
|
||||
struct generic_quantile_finder
|
||||
{
|
||||
typedef typename Dist::value_type value_type;
|
||||
typedef typename Dist::policy_type policy_type;
|
||||
|
||||
generic_quantile_finder(const Dist& d, value_type t, bool c)
|
||||
: dist(d), target(t), comp(c) {}
|
||||
|
||||
value_type operator()(const value_type& x)
|
||||
{
|
||||
return comp ?
|
||||
target - cdf(complement(dist, x))
|
||||
: cdf(dist, x) - target;
|
||||
}
|
||||
|
||||
private:
|
||||
Dist dist;
|
||||
value_type target;
|
||||
bool comp;
|
||||
};
|
||||
|
||||
template <class T, class Policy>
|
||||
inline T check_range_result(const T& x, const Policy& pol, const char* function)
|
||||
{
|
||||
if((x >= 0) && (x < tools::min_value<T>()))
|
||||
return policies::raise_underflow_error<T>(function, 0, pol);
|
||||
if(x <= -tools::max_value<T>())
|
||||
return -policies::raise_overflow_error<T>(function, 0, pol);
|
||||
if(x >= tools::max_value<T>())
|
||||
return policies::raise_overflow_error<T>(function, 0, pol);
|
||||
return x;
|
||||
}
|
||||
|
||||
template <class Dist>
|
||||
typename Dist::value_type generic_quantile(const Dist& dist, const typename Dist::value_type& p, const typename Dist::value_type& guess, bool comp, const char* function)
|
||||
{
|
||||
typedef typename Dist::value_type value_type;
|
||||
typedef typename Dist::policy_type policy_type;
|
||||
typedef typename policies::normalise<
|
||||
policy_type,
|
||||
policies::promote_float<false>,
|
||||
policies::promote_double<false>,
|
||||
policies::discrete_quantile<>,
|
||||
policies::assert_undefined<> >::type forwarding_policy;
|
||||
|
||||
//
|
||||
// Special cases first:
|
||||
//
|
||||
if(p == 0)
|
||||
{
|
||||
return comp
|
||||
? check_range_result(range(dist).second, forwarding_policy(), function)
|
||||
: check_range_result(range(dist).first, forwarding_policy(), function);
|
||||
}
|
||||
if(p == 1)
|
||||
{
|
||||
return !comp
|
||||
? check_range_result(range(dist).second, forwarding_policy(), function)
|
||||
: check_range_result(range(dist).first, forwarding_policy(), function);
|
||||
}
|
||||
|
||||
generic_quantile_finder<Dist> f(dist, p, comp);
|
||||
tools::eps_tolerance<value_type> tol(policies::digits<value_type, forwarding_policy>() - 3);
|
||||
boost::uintmax_t max_iter = policies::get_max_root_iterations<forwarding_policy>();
|
||||
std::pair<value_type, value_type> ir = tools::bracket_and_solve_root(
|
||||
f, guess, value_type(2), true, tol, max_iter, forwarding_policy());
|
||||
value_type result = ir.first + (ir.second - ir.first) / 2;
|
||||
if(max_iter >= policies::get_max_root_iterations<forwarding_policy>())
|
||||
{
|
||||
policies::raise_evaluation_error<value_type>(function, "Unable to locate solution in a reasonable time:"
|
||||
" either there is no answer to quantile"
|
||||
" or the answer is infinite. Current best guess is %1%", result, forwarding_policy());
|
||||
}
|
||||
return result;
|
||||
}
|
||||
|
||||
}}} // namespaces
|
||||
|
||||
#endif // BOOST_MATH_DISTIBUTIONS_DETAIL_GENERIC_QUANTILE_HPP
|
||||
|
||||
@@ -0,0 +1,481 @@
|
||||
// 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, value_type q)
|
||||
: dist(d), target(p < q ? p : q), comp(p < q ? false : true) {}
|
||||
|
||||
value_type operator()(value_type const& x)
|
||||
{
|
||||
return comp ? target - cdf(complement(dist, x)) : 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,
|
||||
const typename Dist::value_type& q,
|
||||
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, q);
|
||||
//
|
||||
// Max bounds of the distribution:
|
||||
//
|
||||
value_type min_bound, max_bound;
|
||||
std::tr1::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;
|
||||
}
|
||||
else
|
||||
{
|
||||
b = a;
|
||||
a = (std::max)(b - 1, value_type(0));
|
||||
if(a < min_bound)
|
||||
a = min_bound;
|
||||
fa = f(a);
|
||||
--count;
|
||||
if(fa == 0)
|
||||
return a;
|
||||
}
|
||||
}
|
||||
}
|
||||
//
|
||||
// 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)(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)(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))
|
||||
{
|
||||
if(count == 0)
|
||||
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))
|
||||
{
|
||||
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)
|
||||
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;
|
||||
//
|
||||
// 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;
|
||||
}
|
||||
//
|
||||
// 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,
|
||||
const typename Dist::value_type& p,
|
||||
const typename Dist::value_type& q,
|
||||
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))
|
||||
return 0;
|
||||
return do_inverse_discrete_quantile(
|
||||
dist,
|
||||
p,
|
||||
q,
|
||||
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,
|
||||
const typename Dist::value_type& q,
|
||||
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
|
||||
if(p <= 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(
|
||||
dist,
|
||||
p,
|
||||
q,
|
||||
(guess < 1 ? value_type(1) : floor(guess)),
|
||||
multiplier,
|
||||
adder,
|
||||
tools::equal_floor(),
|
||||
max_iter));
|
||||
// else:
|
||||
return ceil(do_inverse_discrete_quantile(
|
||||
dist,
|
||||
p,
|
||||
q,
|
||||
ceil(guess),
|
||||
multiplier,
|
||||
adder,
|
||||
tools::equal_ceil(),
|
||||
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& q,
|
||||
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
|
||||
if(p <= 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(
|
||||
dist,
|
||||
p,
|
||||
q,
|
||||
ceil(guess),
|
||||
multiplier,
|
||||
adder,
|
||||
tools::equal_ceil(),
|
||||
max_iter));
|
||||
// else:
|
||||
return floor(do_inverse_discrete_quantile(
|
||||
dist,
|
||||
p,
|
||||
q,
|
||||
(guess < 1 ? value_type(1) : floor(guess)),
|
||||
multiplier,
|
||||
adder,
|
||||
tools::equal_floor(),
|
||||
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& q,
|
||||
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
|
||||
if(p <= pdf(dist, 0))
|
||||
return 0;
|
||||
return floor(do_inverse_discrete_quantile(
|
||||
dist,
|
||||
p,
|
||||
q,
|
||||
(guess < 1 ? value_type(1) : floor(guess)),
|
||||
multiplier,
|
||||
adder,
|
||||
tools::equal_floor(),
|
||||
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& q,
|
||||
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
|
||||
if(p <= pdf(dist, 0))
|
||||
return 0;
|
||||
return ceil(do_inverse_discrete_quantile(
|
||||
dist,
|
||||
p,
|
||||
q,
|
||||
ceil(guess),
|
||||
multiplier,
|
||||
adder,
|
||||
tools::equal_ceil(),
|
||||
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& q,
|
||||
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
|
||||
if(p <= 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(
|
||||
dist,
|
||||
p,
|
||||
q,
|
||||
(guess < 0.5f ? value_type(1.5f) : floor(guess + 0.5f) + 0.5f),
|
||||
multiplier,
|
||||
adder,
|
||||
tools::equal_nearest_integer(),
|
||||
max_iter) + 0.5f);
|
||||
}
|
||||
|
||||
}}} // namespaces
|
||||
|
||||
#endif // BOOST_MATH_DISTRIBUTIONS_DETAIL_INV_DISCRETE_QUANTILE
|
||||
|
||||
|
||||
Reference in New Issue
Block a user