distributions/detail/generic_quantile.hpp

0001 //  Copyright John Maddock 2008.
0002 //  Use, modification and distribution are subject to the
0003 //  Boost Software License, Version 1.0. (See accompanying file
0004 //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
0005
0006 #ifndef BOOST_MATH_DISTIBUTIONS_DETAIL_GENERIC_QUANTILE_HPP
0007 #define BOOST_MATH_DISTIBUTIONS_DETAIL_GENERIC_QUANTILE_HPP
0008
0009 #include <boost/math/tools/config.hpp>
0010 #include <boost/math/tools/tuple.hpp>
0011 #include <boost/math/tools/cstdint.hpp>
0012
0013 namespace boost{ namespace math{ namespace detail{
0014
0015 template <class Dist>
0016 struct generic_quantile_finder
0017 {
0018    using value_type = typename Dist::value_type;
0019    using policy_type = typename Dist::policy_type;
0020
0021    BOOST_MATH_GPU_ENABLED generic_quantile_finder(const Dist& d, value_type t, bool c)
0022       : dist(d), target(t), comp(c) {}
0023
0024    BOOST_MATH_GPU_ENABLED value_type operator()(const value_type& x)
0025    {
0026       return comp ?
0027          value_type(target - cdf(complement(dist, x)))
0028          : value_type(cdf(dist, x) - target);
0029    }
0030
0031 private:
0032    Dist dist;
0033    value_type target;
0034    bool comp;
0035 };
0036
0037 template <class T, class Policy>
0038 BOOST_MATH_GPU_ENABLED inline T check_range_result(const T& x, const Policy& pol, const char* function)
0039 {
0040    if((x >= 0) && (x < tools::min_value<T>()))
0041    {
0042       return policies::raise_underflow_error<T>(function, nullptr, pol);
0043    }
0044    if(x <= -tools::max_value<T>())
0045    {
0046       return -policies::raise_overflow_error<T>(function, nullptr, pol);
0047    }
0048    if(x >= tools::max_value<T>())
0049    {
0050       return policies::raise_overflow_error<T>(function, nullptr, pol);
0051    }
0052    return x;
0053 }
0054
0055 template <class Dist>
0056 BOOST_MATH_GPU_ENABLED 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)
0057 {
0058    using value_type = typename Dist::value_type;
0059    using policy_type = typename Dist::policy_type;
0060    using forwarding_policy = typename policies::normalise<
0061                                                             policy_type,
0062                                                             policies::promote_float<false>,
0063                                                             policies::promote_double<false>,
0064                                                             policies::discrete_quantile<>,
0065                                                             policies::assert_undefined<> >::type;
0066
0067    //
0068    // Special cases first:
0069    //
0070    if(p == 0)
0071    {
0072       return comp
0073       ? check_range_result(range(dist).second, forwarding_policy(), function)
0074       : check_range_result(range(dist).first, forwarding_policy(), function);
0075    }
0076    if(p == 1)
0077    {
0078       return !comp
0079       ? check_range_result(range(dist).second, forwarding_policy(), function)
0080       : check_range_result(range(dist).first, forwarding_policy(), function);
0081    }
0082
0083    generic_quantile_finder<Dist> f(dist, p, comp);
0084    tools::eps_tolerance<value_type> tol(policies::digits<value_type, forwarding_policy>() - 3);
0085    boost::math::uintmax_t max_iter = policies::get_max_root_iterations<forwarding_policy>();
0086    boost::math::pair<value_type, value_type> ir = tools::bracket_and_solve_root(
0087       f, guess, value_type(2), true, tol, max_iter, forwarding_policy());
0088    value_type result = ir.first + (ir.second - ir.first) / 2;
0089    if(max_iter >= policies::get_max_root_iterations<forwarding_policy>())
0090    {
0091       return policies::raise_evaluation_error<value_type>(function, "Unable to locate solution in a reasonable time:" // LCOV_EXCL_LINE
0092          " either there is no answer to quantile or the answer is infinite.  Current best guess is %1%", result, forwarding_policy());  // LCOV_EXCL_LINE
0093    }
0094    return result;
0095 }
0096
0097 }}} // namespaces
0098
0099 #endif // BOOST_MATH_DISTIBUTIONS_DETAIL_GENERIC_QUANTILE_HPP
0100