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