generic_quantile.hpp 3.7 KB

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