annotate DEPENDENCIES/generic/include/boost/math/distributions/detail/generic_mode.hpp @ 133:4acb5d8d80b6 tip

Don't fail environmental check if README.md exists (but .txt and no-suffix don't)
author Chris Cannam
date Tue, 30 Jul 2019 12:25:44 +0100
parents c530137014c0
children
rev   line source
Chris@16 1 // Copyright John Maddock 2008.
Chris@16 2
Chris@16 3 // Use, modification and distribution are subject to the
Chris@16 4 // Boost Software License, Version 1.0.
Chris@16 5 // (See accompanying file LICENSE_1_0.txt
Chris@16 6 // or copy at http://www.boost.org/LICENSE_1_0.txt)
Chris@16 7
Chris@16 8 #ifndef BOOST_MATH_DISTRIBUTIONS_DETAIL_MODE_HPP
Chris@16 9 #define BOOST_MATH_DISTRIBUTIONS_DETAIL_MODE_HPP
Chris@16 10
Chris@16 11 #include <boost/math/tools/minima.hpp> // function minimization for mode
Chris@16 12 #include <boost/math/policies/error_handling.hpp>
Chris@16 13 #include <boost/math/distributions/fwd.hpp>
Chris@16 14
Chris@16 15 namespace boost{ namespace math{ namespace detail{
Chris@16 16
Chris@16 17 template <class Dist>
Chris@16 18 struct pdf_minimizer
Chris@16 19 {
Chris@16 20 pdf_minimizer(const Dist& d)
Chris@16 21 : dist(d) {}
Chris@16 22
Chris@16 23 typename Dist::value_type operator()(const typename Dist::value_type& x)
Chris@16 24 {
Chris@16 25 return -pdf(dist, x);
Chris@16 26 }
Chris@16 27 private:
Chris@16 28 Dist dist;
Chris@16 29 };
Chris@16 30
Chris@16 31 template <class Dist>
Chris@16 32 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 33 {
Chris@16 34 BOOST_MATH_STD_USING
Chris@16 35 typedef typename Dist::value_type value_type;
Chris@16 36 typedef typename Dist::policy_type policy_type;
Chris@16 37 //
Chris@16 38 // Need to begin by bracketing the maxima of the PDF:
Chris@16 39 //
Chris@16 40 value_type maxval;
Chris@16 41 value_type upper_bound = guess;
Chris@16 42 value_type lower_bound;
Chris@16 43 value_type v = pdf(dist, guess);
Chris@16 44 if(v == 0)
Chris@16 45 {
Chris@16 46 //
Chris@16 47 // Oops we don't know how to handle this, or even in which
Chris@16 48 // direction we should move in, treat as an evaluation error:
Chris@16 49 //
Chris@101 50 return policies::raise_evaluation_error(
Chris@16 51 function,
Chris@16 52 "Could not locate a starting location for the search for the mode, original guess was %1%", guess, policy_type());
Chris@16 53 }
Chris@16 54 do
Chris@16 55 {
Chris@16 56 maxval = v;
Chris@16 57 if(step != 0)
Chris@16 58 upper_bound += step;
Chris@16 59 else
Chris@16 60 upper_bound *= 2;
Chris@16 61 v = pdf(dist, upper_bound);
Chris@16 62 }while(maxval < v);
Chris@16 63
Chris@16 64 lower_bound = upper_bound;
Chris@16 65 do
Chris@16 66 {
Chris@16 67 maxval = v;
Chris@16 68 if(step != 0)
Chris@16 69 lower_bound -= step;
Chris@16 70 else
Chris@16 71 lower_bound /= 2;
Chris@16 72 v = pdf(dist, lower_bound);
Chris@16 73 }while(maxval < v);
Chris@16 74
Chris@16 75 boost::uintmax_t max_iter = policies::get_max_root_iterations<policy_type>();
Chris@16 76
Chris@16 77 value_type result = tools::brent_find_minima(
Chris@16 78 pdf_minimizer<Dist>(dist),
Chris@16 79 lower_bound,
Chris@16 80 upper_bound,
Chris@16 81 policies::digits<value_type, policy_type>(),
Chris@16 82 max_iter).first;
Chris@16 83 if(max_iter >= policies::get_max_root_iterations<policy_type>())
Chris@16 84 {
Chris@16 85 return policies::raise_evaluation_error<value_type>(
Chris@16 86 function,
Chris@16 87 "Unable to locate solution in a reasonable time:"
Chris@16 88 " either there is no answer to the mode of the distribution"
Chris@16 89 " or the answer is infinite. Current best guess is %1%", result, policy_type());
Chris@16 90 }
Chris@16 91 return result;
Chris@16 92 }
Chris@16 93 //
Chris@16 94 // As above,but confined to the interval [0,1]:
Chris@16 95 //
Chris@16 96 template <class Dist>
Chris@16 97 typename Dist::value_type generic_find_mode_01(const Dist& dist, typename Dist::value_type guess, const char* function)
Chris@16 98 {
Chris@16 99 BOOST_MATH_STD_USING
Chris@16 100 typedef typename Dist::value_type value_type;
Chris@16 101 typedef typename Dist::policy_type policy_type;
Chris@16 102 //
Chris@16 103 // Need to begin by bracketing the maxima of the PDF:
Chris@16 104 //
Chris@16 105 value_type maxval;
Chris@16 106 value_type upper_bound = guess;
Chris@16 107 value_type lower_bound;
Chris@16 108 value_type v = pdf(dist, guess);
Chris@16 109 do
Chris@16 110 {
Chris@16 111 maxval = v;
Chris@16 112 upper_bound = 1 - (1 - upper_bound) / 2;
Chris@16 113 if(upper_bound == 1)
Chris@16 114 return 1;
Chris@16 115 v = pdf(dist, upper_bound);
Chris@16 116 }while(maxval < v);
Chris@16 117
Chris@16 118 lower_bound = upper_bound;
Chris@16 119 do
Chris@16 120 {
Chris@16 121 maxval = v;
Chris@16 122 lower_bound /= 2;
Chris@16 123 if(lower_bound < tools::min_value<value_type>())
Chris@16 124 return 0;
Chris@16 125 v = pdf(dist, lower_bound);
Chris@16 126 }while(maxval < v);
Chris@16 127
Chris@16 128 boost::uintmax_t max_iter = policies::get_max_root_iterations<policy_type>();
Chris@16 129
Chris@16 130 value_type result = tools::brent_find_minima(
Chris@16 131 pdf_minimizer<Dist>(dist),
Chris@16 132 lower_bound,
Chris@16 133 upper_bound,
Chris@16 134 policies::digits<value_type, policy_type>(),
Chris@16 135 max_iter).first;
Chris@16 136 if(max_iter >= policies::get_max_root_iterations<policy_type>())
Chris@16 137 {
Chris@16 138 return policies::raise_evaluation_error<value_type>(
Chris@16 139 function,
Chris@16 140 "Unable to locate solution in a reasonable time:"
Chris@16 141 " either there is no answer to the mode of the distribution"
Chris@16 142 " or the answer is infinite. Current best guess is %1%", result, policy_type());
Chris@16 143 }
Chris@16 144 return result;
Chris@16 145 }
Chris@16 146
Chris@16 147 }}} // namespaces
Chris@16 148
Chris@16 149 #endif // BOOST_MATH_DISTRIBUTIONS_DETAIL_MODE_HPP