Chris@16
|
1 // Copyright Paul A. Bristow 2007.
|
Chris@16
|
2 // Use, modification and distribution are subject to the
|
Chris@16
|
3 // Boost Software License, Version 1.0. (See accompanying file
|
Chris@16
|
4 // LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
|
Chris@16
|
5
|
Chris@16
|
6 #ifndef BOOST_STATS_rayleigh_HPP
|
Chris@16
|
7 #define BOOST_STATS_rayleigh_HPP
|
Chris@16
|
8
|
Chris@16
|
9 #include <boost/math/distributions/fwd.hpp>
|
Chris@16
|
10 #include <boost/math/constants/constants.hpp>
|
Chris@16
|
11 #include <boost/math/special_functions/log1p.hpp>
|
Chris@16
|
12 #include <boost/math/special_functions/expm1.hpp>
|
Chris@16
|
13 #include <boost/math/distributions/complement.hpp>
|
Chris@16
|
14 #include <boost/math/distributions/detail/common_error_handling.hpp>
|
Chris@16
|
15 #include <boost/config/no_tr1/cmath.hpp>
|
Chris@16
|
16
|
Chris@16
|
17 #ifdef BOOST_MSVC
|
Chris@16
|
18 # pragma warning(push)
|
Chris@16
|
19 # pragma warning(disable: 4702) // unreachable code (return after domain_error throw).
|
Chris@16
|
20 #endif
|
Chris@16
|
21
|
Chris@16
|
22 #include <utility>
|
Chris@16
|
23
|
Chris@16
|
24 namespace boost{ namespace math{
|
Chris@16
|
25
|
Chris@16
|
26 namespace detail
|
Chris@16
|
27 { // Error checks:
|
Chris@16
|
28 template <class RealType, class Policy>
|
Chris@16
|
29 inline bool verify_sigma(const char* function, RealType sigma, RealType* presult, const Policy& pol)
|
Chris@16
|
30 {
|
Chris@16
|
31 if((sigma <= 0) || (!(boost::math::isfinite)(sigma)))
|
Chris@16
|
32 {
|
Chris@16
|
33 *presult = policies::raise_domain_error<RealType>(
|
Chris@16
|
34 function,
|
Chris@16
|
35 "The scale parameter \"sigma\" must be > 0 and finite, but was: %1%.", sigma, pol);
|
Chris@16
|
36 return false;
|
Chris@16
|
37 }
|
Chris@16
|
38 return true;
|
Chris@16
|
39 } // bool verify_sigma
|
Chris@16
|
40
|
Chris@16
|
41 template <class RealType, class Policy>
|
Chris@16
|
42 inline bool verify_rayleigh_x(const char* function, RealType x, RealType* presult, const Policy& pol)
|
Chris@16
|
43 {
|
Chris@16
|
44 if((x < 0) || (boost::math::isnan)(x))
|
Chris@16
|
45 {
|
Chris@16
|
46 *presult = policies::raise_domain_error<RealType>(
|
Chris@16
|
47 function,
|
Chris@16
|
48 "The random variable must be >= 0, but was: %1%.", x, pol);
|
Chris@16
|
49 return false;
|
Chris@16
|
50 }
|
Chris@16
|
51 return true;
|
Chris@16
|
52 } // bool verify_rayleigh_x
|
Chris@16
|
53 } // namespace detail
|
Chris@16
|
54
|
Chris@16
|
55 template <class RealType = double, class Policy = policies::policy<> >
|
Chris@16
|
56 class rayleigh_distribution
|
Chris@16
|
57 {
|
Chris@16
|
58 public:
|
Chris@16
|
59 typedef RealType value_type;
|
Chris@16
|
60 typedef Policy policy_type;
|
Chris@16
|
61
|
Chris@16
|
62 rayleigh_distribution(RealType l_sigma = 1)
|
Chris@16
|
63 : m_sigma(l_sigma)
|
Chris@16
|
64 {
|
Chris@16
|
65 RealType err;
|
Chris@16
|
66 detail::verify_sigma("boost::math::rayleigh_distribution<%1%>::rayleigh_distribution", l_sigma, &err, Policy());
|
Chris@16
|
67 } // rayleigh_distribution
|
Chris@16
|
68
|
Chris@16
|
69 RealType sigma()const
|
Chris@16
|
70 { // Accessor.
|
Chris@16
|
71 return m_sigma;
|
Chris@16
|
72 }
|
Chris@16
|
73
|
Chris@16
|
74 private:
|
Chris@16
|
75 RealType m_sigma;
|
Chris@16
|
76 }; // class rayleigh_distribution
|
Chris@16
|
77
|
Chris@16
|
78 typedef rayleigh_distribution<double> rayleigh;
|
Chris@16
|
79
|
Chris@16
|
80 template <class RealType, class Policy>
|
Chris@16
|
81 inline const std::pair<RealType, RealType> range(const rayleigh_distribution<RealType, Policy>& /*dist*/)
|
Chris@16
|
82 { // Range of permissible values for random variable x.
|
Chris@16
|
83 using boost::math::tools::max_value;
|
Chris@16
|
84 return std::pair<RealType, RealType>(static_cast<RealType>(0), std::numeric_limits<RealType>::has_infinity ? std::numeric_limits<RealType>::infinity() : max_value<RealType>());
|
Chris@16
|
85 }
|
Chris@16
|
86
|
Chris@16
|
87 template <class RealType, class Policy>
|
Chris@16
|
88 inline const std::pair<RealType, RealType> support(const rayleigh_distribution<RealType, Policy>& /*dist*/)
|
Chris@16
|
89 { // Range of supported values for random variable x.
|
Chris@16
|
90 // This is range where cdf rises from 0 to 1, and outside it, the pdf is zero.
|
Chris@16
|
91 using boost::math::tools::max_value;
|
Chris@16
|
92 return std::pair<RealType, RealType>(static_cast<RealType>(0), max_value<RealType>());
|
Chris@16
|
93 }
|
Chris@16
|
94
|
Chris@16
|
95 template <class RealType, class Policy>
|
Chris@16
|
96 inline RealType pdf(const rayleigh_distribution<RealType, Policy>& dist, const RealType& x)
|
Chris@16
|
97 {
|
Chris@16
|
98 BOOST_MATH_STD_USING // for ADL of std function exp.
|
Chris@16
|
99
|
Chris@16
|
100 RealType sigma = dist.sigma();
|
Chris@16
|
101 RealType result = 0;
|
Chris@16
|
102 static const char* function = "boost::math::pdf(const rayleigh_distribution<%1%>&, %1%)";
|
Chris@16
|
103 if(false == detail::verify_sigma(function, sigma, &result, Policy()))
|
Chris@16
|
104 {
|
Chris@16
|
105 return result;
|
Chris@16
|
106 }
|
Chris@16
|
107 if(false == detail::verify_rayleigh_x(function, x, &result, Policy()))
|
Chris@16
|
108 {
|
Chris@16
|
109 return result;
|
Chris@16
|
110 }
|
Chris@16
|
111 if((boost::math::isinf)(x))
|
Chris@16
|
112 {
|
Chris@16
|
113 return 0;
|
Chris@16
|
114 }
|
Chris@16
|
115 RealType sigmasqr = sigma * sigma;
|
Chris@16
|
116 result = x * (exp(-(x * x) / ( 2 * sigmasqr))) / sigmasqr;
|
Chris@16
|
117 return result;
|
Chris@16
|
118 } // pdf
|
Chris@16
|
119
|
Chris@16
|
120 template <class RealType, class Policy>
|
Chris@16
|
121 inline RealType cdf(const rayleigh_distribution<RealType, Policy>& dist, const RealType& x)
|
Chris@16
|
122 {
|
Chris@16
|
123 BOOST_MATH_STD_USING // for ADL of std functions
|
Chris@16
|
124
|
Chris@16
|
125 RealType result = 0;
|
Chris@16
|
126 RealType sigma = dist.sigma();
|
Chris@16
|
127 static const char* function = "boost::math::cdf(const rayleigh_distribution<%1%>&, %1%)";
|
Chris@16
|
128 if(false == detail::verify_sigma(function, sigma, &result, Policy()))
|
Chris@16
|
129 {
|
Chris@16
|
130 return result;
|
Chris@16
|
131 }
|
Chris@16
|
132 if(false == detail::verify_rayleigh_x(function, x, &result, Policy()))
|
Chris@16
|
133 {
|
Chris@16
|
134 return result;
|
Chris@16
|
135 }
|
Chris@16
|
136 result = -boost::math::expm1(-x * x / ( 2 * sigma * sigma), Policy());
|
Chris@16
|
137 return result;
|
Chris@16
|
138 } // cdf
|
Chris@16
|
139
|
Chris@16
|
140 template <class RealType, class Policy>
|
Chris@16
|
141 inline RealType quantile(const rayleigh_distribution<RealType, Policy>& dist, const RealType& p)
|
Chris@16
|
142 {
|
Chris@16
|
143 BOOST_MATH_STD_USING // for ADL of std functions
|
Chris@16
|
144
|
Chris@16
|
145 RealType result = 0;
|
Chris@16
|
146 RealType sigma = dist.sigma();
|
Chris@16
|
147 static const char* function = "boost::math::quantile(const rayleigh_distribution<%1%>&, %1%)";
|
Chris@16
|
148 if(false == detail::verify_sigma(function, sigma, &result, Policy()))
|
Chris@16
|
149 return result;
|
Chris@16
|
150 if(false == detail::check_probability(function, p, &result, Policy()))
|
Chris@16
|
151 return result;
|
Chris@16
|
152
|
Chris@16
|
153 if(p == 0)
|
Chris@16
|
154 {
|
Chris@16
|
155 return 0;
|
Chris@16
|
156 }
|
Chris@16
|
157 if(p == 1)
|
Chris@16
|
158 {
|
Chris@16
|
159 return policies::raise_overflow_error<RealType>(function, 0, Policy());
|
Chris@16
|
160 }
|
Chris@16
|
161 result = sqrt(-2 * sigma * sigma * boost::math::log1p(-p, Policy()));
|
Chris@16
|
162 return result;
|
Chris@16
|
163 } // quantile
|
Chris@16
|
164
|
Chris@16
|
165 template <class RealType, class Policy>
|
Chris@16
|
166 inline RealType cdf(const complemented2_type<rayleigh_distribution<RealType, Policy>, RealType>& c)
|
Chris@16
|
167 {
|
Chris@16
|
168 BOOST_MATH_STD_USING // for ADL of std functions
|
Chris@16
|
169
|
Chris@16
|
170 RealType result = 0;
|
Chris@16
|
171 RealType sigma = c.dist.sigma();
|
Chris@16
|
172 static const char* function = "boost::math::cdf(const rayleigh_distribution<%1%>&, %1%)";
|
Chris@16
|
173 if(false == detail::verify_sigma(function, sigma, &result, Policy()))
|
Chris@16
|
174 {
|
Chris@16
|
175 return result;
|
Chris@16
|
176 }
|
Chris@16
|
177 RealType x = c.param;
|
Chris@16
|
178 if(false == detail::verify_rayleigh_x(function, x, &result, Policy()))
|
Chris@16
|
179 {
|
Chris@16
|
180 return result;
|
Chris@16
|
181 }
|
Chris@101
|
182 RealType ea = x * x / (2 * sigma * sigma);
|
Chris@101
|
183 // Fix for VC11/12 x64 bug in exp(float):
|
Chris@101
|
184 if (ea >= tools::max_value<RealType>())
|
Chris@101
|
185 return 0;
|
Chris@101
|
186 result = exp(-ea);
|
Chris@16
|
187 return result;
|
Chris@16
|
188 } // cdf complement
|
Chris@16
|
189
|
Chris@16
|
190 template <class RealType, class Policy>
|
Chris@16
|
191 inline RealType quantile(const complemented2_type<rayleigh_distribution<RealType, Policy>, RealType>& c)
|
Chris@16
|
192 {
|
Chris@16
|
193 BOOST_MATH_STD_USING // for ADL of std functions, log & sqrt.
|
Chris@16
|
194
|
Chris@16
|
195 RealType result = 0;
|
Chris@16
|
196 RealType sigma = c.dist.sigma();
|
Chris@16
|
197 static const char* function = "boost::math::quantile(const rayleigh_distribution<%1%>&, %1%)";
|
Chris@16
|
198 if(false == detail::verify_sigma(function, sigma, &result, Policy()))
|
Chris@16
|
199 {
|
Chris@16
|
200 return result;
|
Chris@16
|
201 }
|
Chris@16
|
202 RealType q = c.param;
|
Chris@16
|
203 if(false == detail::check_probability(function, q, &result, Policy()))
|
Chris@16
|
204 {
|
Chris@16
|
205 return result;
|
Chris@16
|
206 }
|
Chris@16
|
207 if(q == 1)
|
Chris@16
|
208 {
|
Chris@16
|
209 return 0;
|
Chris@16
|
210 }
|
Chris@16
|
211 if(q == 0)
|
Chris@16
|
212 {
|
Chris@16
|
213 return policies::raise_overflow_error<RealType>(function, 0, Policy());
|
Chris@16
|
214 }
|
Chris@16
|
215 result = sqrt(-2 * sigma * sigma * log(q));
|
Chris@16
|
216 return result;
|
Chris@16
|
217 } // quantile complement
|
Chris@16
|
218
|
Chris@16
|
219 template <class RealType, class Policy>
|
Chris@16
|
220 inline RealType mean(const rayleigh_distribution<RealType, Policy>& dist)
|
Chris@16
|
221 {
|
Chris@16
|
222 RealType result = 0;
|
Chris@16
|
223 RealType sigma = dist.sigma();
|
Chris@16
|
224 static const char* function = "boost::math::mean(const rayleigh_distribution<%1%>&, %1%)";
|
Chris@16
|
225 if(false == detail::verify_sigma(function, sigma, &result, Policy()))
|
Chris@16
|
226 {
|
Chris@16
|
227 return result;
|
Chris@16
|
228 }
|
Chris@16
|
229 using boost::math::constants::root_half_pi;
|
Chris@16
|
230 return sigma * root_half_pi<RealType>();
|
Chris@16
|
231 } // mean
|
Chris@16
|
232
|
Chris@16
|
233 template <class RealType, class Policy>
|
Chris@16
|
234 inline RealType variance(const rayleigh_distribution<RealType, Policy>& dist)
|
Chris@16
|
235 {
|
Chris@16
|
236 RealType result = 0;
|
Chris@16
|
237 RealType sigma = dist.sigma();
|
Chris@16
|
238 static const char* function = "boost::math::variance(const rayleigh_distribution<%1%>&, %1%)";
|
Chris@16
|
239 if(false == detail::verify_sigma(function, sigma, &result, Policy()))
|
Chris@16
|
240 {
|
Chris@16
|
241 return result;
|
Chris@16
|
242 }
|
Chris@16
|
243 using boost::math::constants::four_minus_pi;
|
Chris@16
|
244 return four_minus_pi<RealType>() * sigma * sigma / 2;
|
Chris@16
|
245 } // variance
|
Chris@16
|
246
|
Chris@16
|
247 template <class RealType, class Policy>
|
Chris@16
|
248 inline RealType mode(const rayleigh_distribution<RealType, Policy>& dist)
|
Chris@16
|
249 {
|
Chris@16
|
250 return dist.sigma();
|
Chris@16
|
251 }
|
Chris@16
|
252
|
Chris@16
|
253 template <class RealType, class Policy>
|
Chris@16
|
254 inline RealType median(const rayleigh_distribution<RealType, Policy>& dist)
|
Chris@16
|
255 {
|
Chris@16
|
256 using boost::math::constants::root_ln_four;
|
Chris@16
|
257 return root_ln_four<RealType>() * dist.sigma();
|
Chris@16
|
258 }
|
Chris@16
|
259
|
Chris@16
|
260 template <class RealType, class Policy>
|
Chris@16
|
261 inline RealType skewness(const rayleigh_distribution<RealType, Policy>& /*dist*/)
|
Chris@16
|
262 {
|
Chris@16
|
263 // using namespace boost::math::constants;
|
Chris@16
|
264 return static_cast<RealType>(0.63111065781893713819189935154422777984404221106391L);
|
Chris@16
|
265 // Computed using NTL at 150 bit, about 50 decimal digits.
|
Chris@16
|
266 // return 2 * root_pi<RealType>() * pi_minus_three<RealType>() / pow23_four_minus_pi<RealType>();
|
Chris@16
|
267 }
|
Chris@16
|
268
|
Chris@16
|
269 template <class RealType, class Policy>
|
Chris@16
|
270 inline RealType kurtosis(const rayleigh_distribution<RealType, Policy>& /*dist*/)
|
Chris@16
|
271 {
|
Chris@16
|
272 // using namespace boost::math::constants;
|
Chris@16
|
273 return static_cast<RealType>(3.2450893006876380628486604106197544154170667057995L);
|
Chris@16
|
274 // Computed using NTL at 150 bit, about 50 decimal digits.
|
Chris@16
|
275 // return 3 - (6 * pi<RealType>() * pi<RealType>() - 24 * pi<RealType>() + 16) /
|
Chris@16
|
276 // (four_minus_pi<RealType>() * four_minus_pi<RealType>());
|
Chris@16
|
277 }
|
Chris@16
|
278
|
Chris@16
|
279 template <class RealType, class Policy>
|
Chris@16
|
280 inline RealType kurtosis_excess(const rayleigh_distribution<RealType, Policy>& /*dist*/)
|
Chris@16
|
281 {
|
Chris@16
|
282 //using namespace boost::math::constants;
|
Chris@16
|
283 // Computed using NTL at 150 bit, about 50 decimal digits.
|
Chris@16
|
284 return static_cast<RealType>(0.2450893006876380628486604106197544154170667057995L);
|
Chris@16
|
285 // return -(6 * pi<RealType>() * pi<RealType>() - 24 * pi<RealType>() + 16) /
|
Chris@16
|
286 // (four_minus_pi<RealType>() * four_minus_pi<RealType>());
|
Chris@16
|
287 } // kurtosis
|
Chris@16
|
288
|
Chris@16
|
289 } // namespace math
|
Chris@16
|
290 } // namespace boost
|
Chris@16
|
291
|
Chris@16
|
292 #ifdef BOOST_MSVC
|
Chris@16
|
293 # pragma warning(pop)
|
Chris@16
|
294 #endif
|
Chris@16
|
295
|
Chris@16
|
296 // This include must be at the end, *after* the accessors
|
Chris@16
|
297 // for this distribution have been defined, in order to
|
Chris@16
|
298 // keep compilers that support two-phase lookup happy.
|
Chris@16
|
299 #include <boost/math/distributions/detail/derived_accessors.hpp>
|
Chris@16
|
300
|
Chris@16
|
301 #endif // BOOST_STATS_rayleigh_HPP
|