| // Copyright John Maddock 2008. |
| |
| // Use, modification and distribution are subject to the |
| // Boost Software License, Version 1.0. |
| // (See accompanying file LICENSE_1_0.txt |
| // or copy at http://www.boost.org/LICENSE_1_0.txt) |
| |
| #ifndef BOOST_MATH_DISTRIBUTIONS_DETAIL_MODE_HPP |
| #define BOOST_MATH_DISTRIBUTIONS_DETAIL_MODE_HPP |
| |
| #include <boost/math/tools/minima.hpp> // function minimization for mode |
| #include <boost/math/policies/error_handling.hpp> |
| #include <boost/math/distributions/fwd.hpp> |
| |
| namespace boost{ namespace math{ namespace detail{ |
| |
| template <class Dist> |
| struct pdf_minimizer |
| { |
| pdf_minimizer(const Dist& d) |
| : dist(d) {} |
| |
| typename Dist::value_type operator()(const typename Dist::value_type& x) |
| { |
| return -pdf(dist, x); |
| } |
| private: |
| Dist dist; |
| }; |
| |
| template <class Dist> |
| typename Dist::value_type generic_find_mode(const Dist& dist, typename Dist::value_type guess, const char* function, typename Dist::value_type step = 0) |
| { |
| BOOST_MATH_STD_USING |
| typedef typename Dist::value_type value_type; |
| typedef typename Dist::policy_type policy_type; |
| // |
| // Need to begin by bracketing the maxima of the PDF: |
| // |
| value_type maxval; |
| value_type upper_bound = guess; |
| value_type lower_bound; |
| value_type v = pdf(dist, guess); |
| if(v == 0) |
| { |
| // |
| // Oops we don't know how to handle this, or even in which |
| // direction we should move in, treat as an evaluation error: |
| // |
| policies::raise_evaluation_error( |
| function, |
| "Could not locate a starting location for the search for the mode, original guess was %1%", guess, policy_type()); |
| } |
| do |
| { |
| maxval = v; |
| if(step != 0) |
| upper_bound += step; |
| else |
| upper_bound *= 2; |
| v = pdf(dist, upper_bound); |
| }while(maxval < v); |
| |
| lower_bound = upper_bound; |
| do |
| { |
| maxval = v; |
| if(step != 0) |
| lower_bound -= step; |
| else |
| lower_bound /= 2; |
| v = pdf(dist, lower_bound); |
| }while(maxval < v); |
| |
| boost::uintmax_t max_iter = policies::get_max_root_iterations<policy_type>(); |
| |
| value_type result = tools::brent_find_minima( |
| pdf_minimizer<Dist>(dist), |
| lower_bound, |
| upper_bound, |
| policies::digits<value_type, policy_type>(), |
| max_iter).first; |
| if(max_iter >= policies::get_max_root_iterations<policy_type>()) |
| { |
| return policies::raise_evaluation_error<value_type>( |
| function, |
| "Unable to locate solution in a reasonable time:" |
| " either there is no answer to the mode of the distribution" |
| " or the answer is infinite. Current best guess is %1%", result, policy_type()); |
| } |
| return result; |
| } |
| // |
| // As above,but confined to the interval [0,1]: |
| // |
| template <class Dist> |
| typename Dist::value_type generic_find_mode_01(const Dist& dist, typename Dist::value_type guess, const char* function) |
| { |
| BOOST_MATH_STD_USING |
| typedef typename Dist::value_type value_type; |
| typedef typename Dist::policy_type policy_type; |
| // |
| // Need to begin by bracketing the maxima of the PDF: |
| // |
| value_type maxval; |
| value_type upper_bound = guess; |
| value_type lower_bound; |
| value_type v = pdf(dist, guess); |
| do |
| { |
| maxval = v; |
| upper_bound = 1 - (1 - upper_bound) / 2; |
| if(upper_bound == 1) |
| return 1; |
| v = pdf(dist, upper_bound); |
| }while(maxval < v); |
| |
| lower_bound = upper_bound; |
| do |
| { |
| maxval = v; |
| lower_bound /= 2; |
| if(lower_bound < tools::min_value<value_type>()) |
| return 0; |
| v = pdf(dist, lower_bound); |
| }while(maxval < v); |
| |
| boost::uintmax_t max_iter = policies::get_max_root_iterations<policy_type>(); |
| |
| value_type result = tools::brent_find_minima( |
| pdf_minimizer<Dist>(dist), |
| lower_bound, |
| upper_bound, |
| policies::digits<value_type, policy_type>(), |
| max_iter).first; |
| if(max_iter >= policies::get_max_root_iterations<policy_type>()) |
| { |
| return policies::raise_evaluation_error<value_type>( |
| function, |
| "Unable to locate solution in a reasonable time:" |
| " either there is no answer to the mode of the distribution" |
| " or the answer is infinite. Current best guess is %1%", result, policy_type()); |
| } |
| return result; |
| } |
| |
| }}} // namespaces |
| |
| #endif // BOOST_MATH_DISTRIBUTIONS_DETAIL_MODE_HPP |