101 lines · plain
1// Copyright John Maddock 2008.2// Use, modification and distribution are subject to the3// Boost Software License, Version 1.0. (See accompanying file4// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)5 6#ifndef BOOST_MATH_DISTIBUTIONS_DETAIL_GENERIC_QUANTILE_HPP7#define BOOST_MATH_DISTIBUTIONS_DETAIL_GENERIC_QUANTILE_HPP8 9#include <boost/math/tools/config.hpp>10#include <boost/math/tools/tuple.hpp>11#include <boost/math/tools/cstdint.hpp>12 13namespace boost{ namespace math{ namespace detail{14 15template <class Dist>16struct generic_quantile_finder17{18 using value_type = typename Dist::value_type;19 using policy_type = typename Dist::policy_type;20 21 BOOST_MATH_GPU_ENABLED generic_quantile_finder(const Dist& d, value_type t, bool c)22 : dist(d), target(t), comp(c) {}23 24 BOOST_MATH_GPU_ENABLED value_type operator()(const value_type& x)25 {26 return comp ?27 value_type(target - cdf(complement(dist, x)))28 : value_type(cdf(dist, x) - target);29 }30 31private:32 Dist dist;33 value_type target;34 bool comp;35};36 37template <class T, class Policy>38BOOST_MATH_GPU_ENABLED inline T check_range_result(const T& x, const Policy& pol, const char* function)39{40 if((x >= 0) && (x < tools::min_value<T>()))41 {42 return policies::raise_underflow_error<T>(function, nullptr, pol);43 }44 if(x <= -tools::max_value<T>())45 {46 return -policies::raise_overflow_error<T>(function, nullptr, pol);47 }48 if(x >= tools::max_value<T>())49 {50 return policies::raise_overflow_error<T>(function, nullptr, pol);51 }52 return x;53}54 55template <class Dist>56BOOST_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)57{58 using value_type = typename Dist::value_type;59 using policy_type = typename Dist::policy_type;60 using forwarding_policy = typename policies::normalise<61 policy_type,62 policies::promote_float<false>,63 policies::promote_double<false>,64 policies::discrete_quantile<>,65 policies::assert_undefined<> >::type;66 67 //68 // Special cases first:69 //70 if(p == 0)71 {72 return comp73 ? check_range_result(range(dist).second, forwarding_policy(), function)74 : check_range_result(range(dist).first, forwarding_policy(), function);75 }76 if(p == 1)77 {78 return !comp79 ? check_range_result(range(dist).second, forwarding_policy(), function)80 : check_range_result(range(dist).first, forwarding_policy(), function);81 }82 83 generic_quantile_finder<Dist> f(dist, p, comp);84 tools::eps_tolerance<value_type> tol(policies::digits<value_type, forwarding_policy>() - 3);85 boost::math::uintmax_t max_iter = policies::get_max_root_iterations<forwarding_policy>();86 boost::math::pair<value_type, value_type> ir = tools::bracket_and_solve_root(87 f, guess, value_type(2), true, tol, max_iter, forwarding_policy());88 value_type result = ir.first + (ir.second - ir.first) / 2;89 if(max_iter >= policies::get_max_root_iterations<forwarding_policy>())90 {91 return policies::raise_evaluation_error<value_type>(function, "Unable to locate solution in a reasonable time:" // LCOV_EXCL_LINE92 " either there is no answer to quantile or the answer is infinite. Current best guess is %1%", result, forwarding_policy()); // LCOV_EXCL_LINE93 }94 return result;95}96 97}}} // namespaces98 99#endif // BOOST_MATH_DISTIBUTIONS_DETAIL_GENERIC_QUANTILE_HPP100 101