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

source code of boost/libs/math/include/boost/math/distributions/detail/generic_mode.hpp