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@16
|
182 result = exp(-x * x / ( 2 * sigma * sigma));
|
Chris@16
|
183 return result;
|
Chris@16
|
184 } // cdf complement
|
Chris@16
|
185
|
Chris@16
|
186 template <class RealType, class Policy>
|
Chris@16
|
187 inline RealType quantile(const complemented2_type<rayleigh_distribution<RealType, Policy>, RealType>& c)
|
Chris@16
|
188 {
|
Chris@16
|
189 BOOST_MATH_STD_USING // for ADL of std functions, log & sqrt.
|
Chris@16
|
190
|
Chris@16
|
191 RealType result = 0;
|
Chris@16
|
192 RealType sigma = c.dist.sigma();
|
Chris@16
|
193 static const char* function = "boost::math::quantile(const rayleigh_distribution<%1%>&, %1%)";
|
Chris@16
|
194 if(false == detail::verify_sigma(function, sigma, &result, Policy()))
|
Chris@16
|
195 {
|
Chris@16
|
196 return result;
|
Chris@16
|
197 }
|
Chris@16
|
198 RealType q = c.param;
|
Chris@16
|
199 if(false == detail::check_probability(function, q, &result, Policy()))
|
Chris@16
|
200 {
|
Chris@16
|
201 return result;
|
Chris@16
|
202 }
|
Chris@16
|
203 if(q == 1)
|
Chris@16
|
204 {
|
Chris@16
|
205 return 0;
|
Chris@16
|
206 }
|
Chris@16
|
207 if(q == 0)
|
Chris@16
|
208 {
|
Chris@16
|
209 return policies::raise_overflow_error<RealType>(function, 0, Policy());
|
Chris@16
|
210 }
|
Chris@16
|
211 result = sqrt(-2 * sigma * sigma * log(q));
|
Chris@16
|
212 return result;
|
Chris@16
|
213 } // quantile complement
|
Chris@16
|
214
|
Chris@16
|
215 template <class RealType, class Policy>
|
Chris@16
|
216 inline RealType mean(const rayleigh_distribution<RealType, Policy>& dist)
|
Chris@16
|
217 {
|
Chris@16
|
218 RealType result = 0;
|
Chris@16
|
219 RealType sigma = dist.sigma();
|
Chris@16
|
220 static const char* function = "boost::math::mean(const rayleigh_distribution<%1%>&, %1%)";
|
Chris@16
|
221 if(false == detail::verify_sigma(function, sigma, &result, Policy()))
|
Chris@16
|
222 {
|
Chris@16
|
223 return result;
|
Chris@16
|
224 }
|
Chris@16
|
225 using boost::math::constants::root_half_pi;
|
Chris@16
|
226 return sigma * root_half_pi<RealType>();
|
Chris@16
|
227 } // mean
|
Chris@16
|
228
|
Chris@16
|
229 template <class RealType, class Policy>
|
Chris@16
|
230 inline RealType variance(const rayleigh_distribution<RealType, Policy>& dist)
|
Chris@16
|
231 {
|
Chris@16
|
232 RealType result = 0;
|
Chris@16
|
233 RealType sigma = dist.sigma();
|
Chris@16
|
234 static const char* function = "boost::math::variance(const rayleigh_distribution<%1%>&, %1%)";
|
Chris@16
|
235 if(false == detail::verify_sigma(function, sigma, &result, Policy()))
|
Chris@16
|
236 {
|
Chris@16
|
237 return result;
|
Chris@16
|
238 }
|
Chris@16
|
239 using boost::math::constants::four_minus_pi;
|
Chris@16
|
240 return four_minus_pi<RealType>() * sigma * sigma / 2;
|
Chris@16
|
241 } // variance
|
Chris@16
|
242
|
Chris@16
|
243 template <class RealType, class Policy>
|
Chris@16
|
244 inline RealType mode(const rayleigh_distribution<RealType, Policy>& dist)
|
Chris@16
|
245 {
|
Chris@16
|
246 return dist.sigma();
|
Chris@16
|
247 }
|
Chris@16
|
248
|
Chris@16
|
249 template <class RealType, class Policy>
|
Chris@16
|
250 inline RealType median(const rayleigh_distribution<RealType, Policy>& dist)
|
Chris@16
|
251 {
|
Chris@16
|
252 using boost::math::constants::root_ln_four;
|
Chris@16
|
253 return root_ln_four<RealType>() * dist.sigma();
|
Chris@16
|
254 }
|
Chris@16
|
255
|
Chris@16
|
256 template <class RealType, class Policy>
|
Chris@16
|
257 inline RealType skewness(const rayleigh_distribution<RealType, Policy>& /*dist*/)
|
Chris@16
|
258 {
|
Chris@16
|
259 // using namespace boost::math::constants;
|
Chris@16
|
260 return static_cast<RealType>(0.63111065781893713819189935154422777984404221106391L);
|
Chris@16
|
261 // Computed using NTL at 150 bit, about 50 decimal digits.
|
Chris@16
|
262 // return 2 * root_pi<RealType>() * pi_minus_three<RealType>() / pow23_four_minus_pi<RealType>();
|
Chris@16
|
263 }
|
Chris@16
|
264
|
Chris@16
|
265 template <class RealType, class Policy>
|
Chris@16
|
266 inline RealType kurtosis(const rayleigh_distribution<RealType, Policy>& /*dist*/)
|
Chris@16
|
267 {
|
Chris@16
|
268 // using namespace boost::math::constants;
|
Chris@16
|
269 return static_cast<RealType>(3.2450893006876380628486604106197544154170667057995L);
|
Chris@16
|
270 // Computed using NTL at 150 bit, about 50 decimal digits.
|
Chris@16
|
271 // return 3 - (6 * pi<RealType>() * pi<RealType>() - 24 * pi<RealType>() + 16) /
|
Chris@16
|
272 // (four_minus_pi<RealType>() * four_minus_pi<RealType>());
|
Chris@16
|
273 }
|
Chris@16
|
274
|
Chris@16
|
275 template <class RealType, class Policy>
|
Chris@16
|
276 inline RealType kurtosis_excess(const rayleigh_distribution<RealType, Policy>& /*dist*/)
|
Chris@16
|
277 {
|
Chris@16
|
278 //using namespace boost::math::constants;
|
Chris@16
|
279 // Computed using NTL at 150 bit, about 50 decimal digits.
|
Chris@16
|
280 return static_cast<RealType>(0.2450893006876380628486604106197544154170667057995L);
|
Chris@16
|
281 // return -(6 * pi<RealType>() * pi<RealType>() - 24 * pi<RealType>() + 16) /
|
Chris@16
|
282 // (four_minus_pi<RealType>() * four_minus_pi<RealType>());
|
Chris@16
|
283 } // kurtosis
|
Chris@16
|
284
|
Chris@16
|
285 } // namespace math
|
Chris@16
|
286 } // namespace boost
|
Chris@16
|
287
|
Chris@16
|
288 #ifdef BOOST_MSVC
|
Chris@16
|
289 # pragma warning(pop)
|
Chris@16
|
290 #endif
|
Chris@16
|
291
|
Chris@16
|
292 // This include must be at the end, *after* the accessors
|
Chris@16
|
293 // for this distribution have been defined, in order to
|
Chris@16
|
294 // keep compilers that support two-phase lookup happy.
|
Chris@16
|
295 #include <boost/math/distributions/detail/derived_accessors.hpp>
|
Chris@16
|
296
|
Chris@16
|
297 #endif // BOOST_STATS_rayleigh_HPP
|