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(
51 function,
52 "Could not locate a starting location for the search for the mode, original guess was %1%", guess, policy_type());
53 }
54 do
55 {
56 maxval = v;
57 if(step != 0)
58 upper_bound += step;
59 else
60 upper_bound *= 2;
61 v = pdf(dist, upper_bound);
62 }while(maxval < v);
63
64 lower_bound = upper_bound;
65 do
66 {
67 maxval = v;
68 if(step != 0)
69 lower_bound -= step;
70 else
71 lower_bound /= 2;
72 v = pdf(dist, lower_bound);
73 }while(maxval < v);
74
75 boost::uintmax_t max_iter = policies::get_max_root_iterations<policy_type>();
76
77 value_type result = tools::brent_find_minima(
78 pdf_minimizer<Dist>(dist),
79 lower_bound,
80 upper_bound,
81 policies::digits<value_type, policy_type>(),
82 max_iter).first;
83 if(max_iter >= policies::get_max_root_iterations<policy_type>())
84 {
85 return policies::raise_evaluation_error<value_type>(
86 function,
87 "Unable to locate solution in a reasonable time:"
88 " either there is no answer to the mode of the distribution"
89 " or the answer is infinite. Current best guess is %1%", result, policy_type());
90 }
91 return result;
92}
93//
94// As above,but confined to the interval [0,1]:
95//
96template <class Dist>
97typename Dist::value_type generic_find_mode_01(const Dist& dist, typename Dist::value_type guess, const char* function)
98{
99 BOOST_MATH_STD_USING
100 typedef typename Dist::value_type value_type;
101 typedef typename Dist::policy_type policy_type;
102 //
103 // Need to begin by bracketing the maxima of the PDF:
104 //
105 value_type maxval;
106 value_type upper_bound = guess;
107 value_type lower_bound;
108 value_type v = pdf(dist, guess);
109 do
110 {
111 maxval = v;
112 upper_bound = 1 - (1 - upper_bound) / 2;
113 if(upper_bound == 1)
114 return 1;
115 v = pdf(dist, upper_bound);
116 }while(maxval < v);
117
118 lower_bound = upper_bound;
119 do
120 {
121 maxval = v;
122 lower_bound /= 2;
123 if(lower_bound < tools::min_value<value_type>())
124 return 0;
125 v = pdf(dist, lower_bound);
126 }while(maxval < v);
127
128 boost::uintmax_t max_iter = policies::get_max_root_iterations<policy_type>();
129
130 value_type result = tools::brent_find_minima(
131 pdf_minimizer<Dist>(dist),
132 lower_bound,
133 upper_bound,
134 policies::digits<value_type, policy_type>(),
135 max_iter).first;
136 if(max_iter >= policies::get_max_root_iterations<policy_type>())
137 {
138 return policies::raise_evaluation_error<value_type>(
139 function,
140 "Unable to locate solution in a reasonable time:"
141 " either there is no answer to the mode of the distribution"
142 " or the answer is infinite. Current best guess is %1%", result, policy_type());
143 }
144 return result;
145}
146
147}}} // namespaces
148
149#endif // BOOST_MATH_DISTRIBUTIONS_DETAIL_MODE_HPP
150

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