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