generic_quantile.hpp 3.1 KB

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