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 | |
15 | namespace boost{ namespace math{ namespace detail{ |
16 | |
17 | template <class Dist> |
18 | struct 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 | } |
27 | private: |
28 | Dist dist; |
29 | }; |
30 | |
31 | template <class Dist> |
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) |
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 | // |
96 | template <class Dist> |
97 | typename 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 | |