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