146 lines · plain
1// Copyright John Maddock 2008.2 3// Use, modification and distribution are subject to the4// Boost Software License, Version 1.0.5// (See accompanying file LICENSE_1_0.txt6// or copy at http://www.boost.org/LICENSE_1_0.txt)7 8#ifndef BOOST_MATH_DISTRIBUTIONS_DETAIL_MODE_HPP9#define BOOST_MATH_DISTRIBUTIONS_DETAIL_MODE_HPP10 11#include <boost/math/tools/config.hpp>12#include <boost/math/tools/cstdint.hpp>13#include <boost/math/tools/minima.hpp> // function minimization for mode14#include <boost/math/policies/error_handling.hpp>15#include <boost/math/distributions/fwd.hpp>16#include <boost/math/policies/policy.hpp>17 18namespace boost{ namespace math{ namespace detail{19 20template <class Dist>21struct pdf_minimizer22{23 BOOST_MATH_GPU_ENABLED pdf_minimizer(const Dist& d)24 : dist(d) {}25 26 BOOST_MATH_GPU_ENABLED typename Dist::value_type operator()(const typename Dist::value_type& x)27 {28 return -pdf(dist, x);29 }30private:31 Dist dist;32};33 34template <class Dist>35BOOST_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)36{37 BOOST_MATH_STD_USING38 typedef typename Dist::value_type value_type;39 typedef typename Dist::policy_type policy_type;40 //41 // Need to begin by bracketing the maxima of the PDF:42 //43 value_type maxval;44 value_type upper_bound = guess;45 value_type lower_bound;46 value_type v = pdf(dist, guess);47 if(v == 0)48 {49 //50 // Oops we don't know how to handle this, or even in which51 // direction we should move in, treat as an evaluation error:52 //53 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_LINE54 }55 do56 {57 maxval = v;58 if(step != 0)59 upper_bound += step;60 else61 upper_bound *= 2;62 v = pdf(dist, upper_bound);63 }while(maxval < v);64 65 lower_bound = upper_bound;66 do67 {68 maxval = v;69 if(step != 0)70 lower_bound -= step;71 else72 lower_bound /= 2;73 v = pdf(dist, lower_bound);74 }while(maxval < v);75 76 boost::math::uintmax_t max_iter = policies::get_max_root_iterations<policy_type>();77 78 value_type result = tools::brent_find_minima(79 pdf_minimizer<Dist>(dist), 80 lower_bound, 81 upper_bound, 82 policies::digits<value_type, policy_type>(), 83 max_iter).first;84 if(max_iter >= policies::get_max_root_iterations<policy_type>())85 {86 return policies::raise_evaluation_error<value_type>(function, // LCOV_EXCL_LINE87 "Unable to locate solution in a reasonable time: either there is no answer to the mode of the distribution" // LCOV_EXCL_LINE88 " or the answer is infinite. Current best guess is %1%", result, policy_type()); // LCOV_EXCL_LINE89 }90 return result;91}92//93// As above,but confined to the interval [0,1]:94//95template <class Dist>96BOOST_MATH_GPU_ENABLED typename Dist::value_type generic_find_mode_01(const Dist& dist, typename Dist::value_type guess, const char* function)97{98 BOOST_MATH_STD_USING99 typedef typename Dist::value_type value_type;100 typedef typename Dist::policy_type policy_type;101 //102 // Need to begin by bracketing the maxima of the PDF:103 //104 value_type maxval;105 value_type upper_bound = guess;106 value_type lower_bound;107 value_type v = pdf(dist, guess);108 do109 {110 maxval = v;111 upper_bound = 1 - (1 - upper_bound) / 2;112 if(upper_bound == 1)113 return 1;114 v = pdf(dist, upper_bound);115 }while(maxval < v);116 117 lower_bound = upper_bound;118 do119 {120 maxval = v;121 lower_bound /= 2;122 if(lower_bound < tools::min_value<value_type>())123 return 0;124 v = pdf(dist, lower_bound);125 }while(maxval < v);126 127 boost::math::uintmax_t max_iter = policies::get_max_root_iterations<policy_type>();128 129 value_type result = tools::brent_find_minima(130 pdf_minimizer<Dist>(dist), 131 lower_bound, 132 upper_bound, 133 policies::digits<value_type, policy_type>(), 134 max_iter).first;135 if(max_iter >= policies::get_max_root_iterations<policy_type>())136 {137 return policies::raise_evaluation_error<value_type>(function, "Unable to locate solution in a reasonable time:" // LCOV_EXCL_LINE138 " 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_LINE139 }140 return result;141}142 143}}} // namespaces144 145#endif // BOOST_MATH_DISTRIBUTIONS_DETAIL_MODE_HPP146