generic_mode.hpp 4.5 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146
  1. // Copyright John Maddock 2008.
  2. // Use, modification and distribution are subject to the
  3. // Boost Software License, Version 1.0.
  4. // (See accompanying file LICENSE_1_0.txt
  5. // or copy at http://www.boost.org/LICENSE_1_0.txt)
  6. #ifndef BOOST_MATH_DISTRIBUTIONS_DETAIL_MODE_HPP
  7. #define BOOST_MATH_DISTRIBUTIONS_DETAIL_MODE_HPP
  8. #include <boost/math/tools/config.hpp>
  9. #include <boost/math/tools/cstdint.hpp>
  10. #include <boost/math/tools/minima.hpp> // function minimization for mode
  11. #include <boost/math/policies/error_handling.hpp>
  12. #include <boost/math/distributions/fwd.hpp>
  13. #include <boost/math/policies/policy.hpp>
  14. namespace boost{ namespace math{ namespace detail{
  15. template <class Dist>
  16. struct pdf_minimizer
  17. {
  18. BOOST_MATH_GPU_ENABLED pdf_minimizer(const Dist& d)
  19. : dist(d) {}
  20. BOOST_MATH_GPU_ENABLED typename Dist::value_type operator()(const typename Dist::value_type& x)
  21. {
  22. return -pdf(dist, x);
  23. }
  24. private:
  25. Dist dist;
  26. };
  27. template <class Dist>
  28. BOOST_MATH_GPU_ENABLED typename Dist::value_type generic_find_mode(const Dist& dist, typename Dist::value_type guess, const char* function, typename Dist::value_type step = 0)
  29. {
  30. BOOST_MATH_STD_USING
  31. typedef typename Dist::value_type value_type;
  32. typedef typename Dist::policy_type policy_type;
  33. //
  34. // Need to begin by bracketing the maxima of the PDF:
  35. //
  36. value_type maxval;
  37. value_type upper_bound = guess;
  38. value_type lower_bound;
  39. value_type v = pdf(dist, guess);
  40. if(v == 0)
  41. {
  42. //
  43. // Oops we don't know how to handle this, or even in which
  44. // direction we should move in, treat as an evaluation error:
  45. //
  46. return policies::raise_evaluation_error(function, "Could not locate a starting location for the search for the mode, original guess was %1%", guess, policy_type()); // LCOV_EXCL_LINE
  47. }
  48. do
  49. {
  50. maxval = v;
  51. if(step != 0)
  52. upper_bound += step;
  53. else
  54. upper_bound *= 2;
  55. v = pdf(dist, upper_bound);
  56. }while(maxval < v);
  57. lower_bound = upper_bound;
  58. do
  59. {
  60. maxval = v;
  61. if(step != 0)
  62. lower_bound -= step;
  63. else
  64. lower_bound /= 2;
  65. v = pdf(dist, lower_bound);
  66. }while(maxval < v);
  67. boost::math::uintmax_t max_iter = policies::get_max_root_iterations<policy_type>();
  68. value_type result = tools::brent_find_minima(
  69. pdf_minimizer<Dist>(dist),
  70. lower_bound,
  71. upper_bound,
  72. policies::digits<value_type, policy_type>(),
  73. max_iter).first;
  74. if(max_iter >= policies::get_max_root_iterations<policy_type>())
  75. {
  76. return policies::raise_evaluation_error<value_type>(function, // LCOV_EXCL_LINE
  77. "Unable to locate solution in a reasonable time: either there is no answer to the mode of the distribution" // LCOV_EXCL_LINE
  78. " or the answer is infinite. Current best guess is %1%", result, policy_type()); // LCOV_EXCL_LINE
  79. }
  80. return result;
  81. }
  82. //
  83. // As above,but confined to the interval [0,1]:
  84. //
  85. template <class Dist>
  86. BOOST_MATH_GPU_ENABLED typename Dist::value_type generic_find_mode_01(const Dist& dist, typename Dist::value_type guess, const char* function)
  87. {
  88. BOOST_MATH_STD_USING
  89. typedef typename Dist::value_type value_type;
  90. typedef typename Dist::policy_type policy_type;
  91. //
  92. // Need to begin by bracketing the maxima of the PDF:
  93. //
  94. value_type maxval;
  95. value_type upper_bound = guess;
  96. value_type lower_bound;
  97. value_type v = pdf(dist, guess);
  98. do
  99. {
  100. maxval = v;
  101. upper_bound = 1 - (1 - upper_bound) / 2;
  102. if(upper_bound == 1)
  103. return 1;
  104. v = pdf(dist, upper_bound);
  105. }while(maxval < v);
  106. lower_bound = upper_bound;
  107. do
  108. {
  109. maxval = v;
  110. lower_bound /= 2;
  111. if(lower_bound < tools::min_value<value_type>())
  112. return 0;
  113. v = pdf(dist, lower_bound);
  114. }while(maxval < v);
  115. boost::math::uintmax_t max_iter = policies::get_max_root_iterations<policy_type>();
  116. value_type result = tools::brent_find_minima(
  117. pdf_minimizer<Dist>(dist),
  118. lower_bound,
  119. upper_bound,
  120. policies::digits<value_type, policy_type>(),
  121. max_iter).first;
  122. if(max_iter >= policies::get_max_root_iterations<policy_type>())
  123. {
  124. return policies::raise_evaluation_error<value_type>(function, "Unable to locate solution in a reasonable time:" // LCOV_EXCL_LINE
  125. " either there is no answer to the mode of the distribution or the answer is infinite. Current best guess is %1%", result, policy_type()); // LCOV_EXCL_LINE
  126. }
  127. return result;
  128. }
  129. }}} // namespaces
  130. #endif // BOOST_MATH_DISTRIBUTIONS_DETAIL_MODE_HPP