distributions/detail/generic_mode.hpp

0001 // Copyright John Maddock 2008.
0002
0003 // Use, modification and distribution are subject to the
0004 // Boost Software License, Version 1.0.
0005 // (See accompanying file LICENSE_1_0.txt
0006 // or copy at http://www.boost.org/LICENSE_1_0.txt)
0007
0008 #ifndef BOOST_MATH_DISTRIBUTIONS_DETAIL_MODE_HPP
0009 #define BOOST_MATH_DISTRIBUTIONS_DETAIL_MODE_HPP
0010
0011 #include <boost/math/tools/config.hpp>
0012 #include <boost/math/tools/cstdint.hpp>
0013 #include <boost/math/tools/minima.hpp> // function minimization for mode
0014 #include <boost/math/policies/error_handling.hpp>
0015 #include <boost/math/distributions/fwd.hpp>
0016 #include <boost/math/policies/policy.hpp>
0017
0018 namespace boost{ namespace math{ namespace detail{
0019
0020 template <class Dist>
0021 struct pdf_minimizer
0022 {
0023    BOOST_MATH_GPU_ENABLED pdf_minimizer(const Dist& d)
0024       : dist(d) {}
0025
0026    BOOST_MATH_GPU_ENABLED typename Dist::value_type operator()(const typename Dist::value_type& x)
0027    {
0028       return -pdf(dist, x);
0029    }
0030 private:
0031    Dist dist;
0032 };
0033
0034 template <class Dist>
0035 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)
0036 {
0037    BOOST_MATH_STD_USING
0038    typedef typename Dist::value_type value_type;
0039    typedef typename Dist::policy_type policy_type;
0040    //
0041    // Need to begin by bracketing the maxima of the PDF:
0042    //
0043    value_type maxval;
0044    value_type upper_bound = guess;
0045    value_type lower_bound;
0046    value_type v = pdf(dist, guess);
0047    if(v == 0)
0048    {
0049       //
0050       // Oops we don't know how to handle this, or even in which
0051       // direction we should move in, treat as an evaluation error:
0052       //
0053       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
0054    }
0055    do
0056    {
0057       maxval = v;
0058       if(step != 0)
0059          upper_bound += step;
0060       else
0061          upper_bound *= 2;
0062       v = pdf(dist, upper_bound);
0063    }while(maxval < v);
0064
0065    lower_bound = upper_bound;
0066    do
0067    {
0068       maxval = v;
0069       if(step != 0)
0070          lower_bound -= step;
0071       else
0072          lower_bound /= 2;
0073       v = pdf(dist, lower_bound);
0074    }while(maxval < v);
0075
0076    boost::math::uintmax_t max_iter = policies::get_max_root_iterations<policy_type>();
0077
0078    value_type result = tools::brent_find_minima(
0079       pdf_minimizer<Dist>(dist),
0080       lower_bound,
0081       upper_bound,
0082       policies::digits<value_type, policy_type>(),
0083       max_iter).first;
0084    if(max_iter >= policies::get_max_root_iterations<policy_type>())
0085    {
0086       return policies::raise_evaluation_error<value_type>(function,   // LCOV_EXCL_LINE
0087          "Unable to locate solution in a reasonable time: either there is no answer to the mode of the distribution"  // LCOV_EXCL_LINE
0088          " or the answer is infinite.  Current best guess is %1%", result, policy_type());  // LCOV_EXCL_LINE
0089    }
0090    return result;
0091 }
0092 //
0093 // As above,but confined to the interval [0,1]:
0094 //
0095 template <class Dist>
0096 BOOST_MATH_GPU_ENABLED typename Dist::value_type generic_find_mode_01(const Dist& dist, typename Dist::value_type guess, const char* function)
0097 {
0098    BOOST_MATH_STD_USING
0099    typedef typename Dist::value_type value_type;
0100    typedef typename Dist::policy_type policy_type;
0101    //
0102    // Need to begin by bracketing the maxima of the PDF:
0103    //
0104    value_type maxval;
0105    value_type upper_bound = guess;
0106    value_type lower_bound;
0107    value_type v = pdf(dist, guess);
0108    do
0109    {
0110       maxval = v;
0111       upper_bound = 1 - (1 - upper_bound) / 2;
0112       if(upper_bound == 1)
0113          return 1;
0114       v = pdf(dist, upper_bound);
0115    }while(maxval < v);
0116
0117    lower_bound = upper_bound;
0118    do
0119    {
0120       maxval = v;
0121       lower_bound /= 2;
0122       if(lower_bound < tools::min_value<value_type>())
0123          return 0;
0124       v = pdf(dist, lower_bound);
0125    }while(maxval < v);
0126
0127    boost::math::uintmax_t max_iter = policies::get_max_root_iterations<policy_type>();
0128
0129    value_type result = tools::brent_find_minima(
0130       pdf_minimizer<Dist>(dist),
0131       lower_bound,
0132       upper_bound,
0133       policies::digits<value_type, policy_type>(),
0134       max_iter).first;
0135    if(max_iter >= policies::get_max_root_iterations<policy_type>())
0136    {
0137       return policies::raise_evaluation_error<value_type>(function, "Unable to locate solution in a reasonable time:" // LCOV_EXCL_LINE
0138          " 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
0139    }
0140    return result;
0141 }
0142
0143 }}} // namespaces
0144
0145 #endif // BOOST_MATH_DISTRIBUTIONS_DETAIL_MODE_HPP