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

source code of include/boost/math/distributions/detail/generic_quantile.hpp