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