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