Mercurial > hg > sv-dependency-builds
comparison any/include/boost/math/distributions/arcsine.hpp @ 160:cff480c41f97
Add some cross-platform Boost headers
author | Chris Cannam <cannam@all-day-breakfast.com> |
---|---|
date | Sat, 16 Feb 2019 16:31:25 +0000 |
parents | |
children |
comparison
equal
deleted
inserted
replaced
159:f4b37539fcc7 | 160:cff480c41f97 |
---|---|
1 // boost/math/distributions/arcsine.hpp | |
2 | |
3 // Copyright John Maddock 2014. | |
4 // Copyright Paul A. Bristow 2014. | |
5 | |
6 // Use, modification and distribution are subject to the | |
7 // Boost Software License, Version 1.0. | |
8 // (See accompanying file LICENSE_1_0.txt | |
9 // or copy at http://www.boost.org/LICENSE_1_0.txt) | |
10 | |
11 // http://en.wikipedia.org/wiki/arcsine_distribution | |
12 | |
13 // The arcsine Distribution is a continuous probability distribution. | |
14 // http://en.wikipedia.org/wiki/Arcsine_distribution | |
15 // http://www.wolframalpha.com/input/?i=ArcSinDistribution | |
16 | |
17 // Standard arcsine distribution is a special case of beta distribution with both a & b = one half, | |
18 // and 0 <= x <= 1. | |
19 | |
20 // It is generalized to include any bounded support a <= x <= b from 0 <= x <= 1 | |
21 // by Wolfram and Wikipedia, | |
22 // but using location and scale parameters by | |
23 // Virtual Laboratories in Probability and Statistics http://www.math.uah.edu/stat/index.html | |
24 // http://www.math.uah.edu/stat/special/Arcsine.html | |
25 // The end-point version is simpler and more obvious, so we implement that. | |
26 // TODO Perhaps provide location and scale functions? | |
27 | |
28 | |
29 #ifndef BOOST_MATH_DIST_ARCSINE_HPP | |
30 #define BOOST_MATH_DIST_ARCSINE_HPP | |
31 | |
32 #include <boost/math/distributions/fwd.hpp> | |
33 #include <boost/math/distributions/complement.hpp> // complements. | |
34 #include <boost/math/distributions/detail/common_error_handling.hpp> // error checks. | |
35 #include <boost/math/constants/constants.hpp> | |
36 | |
37 #include <boost/math/special_functions/fpclassify.hpp> // isnan. | |
38 | |
39 #if defined (BOOST_MSVC) | |
40 # pragma warning(push) | |
41 # pragma warning(disable: 4702) // Unreachable code, | |
42 // in domain_error_imp in error_handling. | |
43 #endif | |
44 | |
45 #include <utility> | |
46 #include <exception> // For std::domain_error. | |
47 | |
48 namespace boost | |
49 { | |
50 namespace math | |
51 { | |
52 namespace arcsine_detail | |
53 { | |
54 // Common error checking routines for arcsine distribution functions: | |
55 // Duplicating for x_min and x_max provides specific error messages. | |
56 template <class RealType, class Policy> | |
57 inline bool check_x_min(const char* function, const RealType& x, RealType* result, const Policy& pol) | |
58 { | |
59 if (!(boost::math::isfinite)(x)) | |
60 { | |
61 *result = policies::raise_domain_error<RealType>( | |
62 function, | |
63 "x_min argument is %1%, but must be finite !", x, pol); | |
64 return false; | |
65 } | |
66 return true; | |
67 } // bool check_x_min | |
68 | |
69 template <class RealType, class Policy> | |
70 inline bool check_x_max(const char* function, const RealType& x, RealType* result, const Policy& pol) | |
71 { | |
72 if (!(boost::math::isfinite)(x)) | |
73 { | |
74 *result = policies::raise_domain_error<RealType>( | |
75 function, | |
76 "x_max argument is %1%, but must be finite !", x, pol); | |
77 return false; | |
78 } | |
79 return true; | |
80 } // bool check_x_max | |
81 | |
82 | |
83 template <class RealType, class Policy> | |
84 inline bool check_x_minmax(const char* function, const RealType& x_min, const RealType& x_max, RealType* result, const Policy& pol) | |
85 { // Check x_min < x_max | |
86 if (x_min >= x_max) | |
87 { | |
88 std::string msg = "x_max argument is %1%, but must be > x_min = " + lexical_cast<std::string>(x_min) + "!"; | |
89 *result = policies::raise_domain_error<RealType>( | |
90 function, | |
91 msg.c_str(), x_max, pol); | |
92 // "x_max argument is %1%, but must be > x_min !", x_max, pol); | |
93 // "x_max argument is %1%, but must be > x_min %2!", x_max, x_min, pol); would be better. | |
94 // But would require replication of all helpers functions in /policies/error_handling.hpp for two values, | |
95 // as well as two value versions of raise_error, raise_domain_error and do_format ... | |
96 // so use slightly hacky lexical_cast to string instead. | |
97 return false; | |
98 } | |
99 return true; | |
100 } // bool check_x_minmax | |
101 | |
102 template <class RealType, class Policy> | |
103 inline bool check_prob(const char* function, const RealType& p, RealType* result, const Policy& pol) | |
104 { | |
105 if ((p < 0) || (p > 1) || !(boost::math::isfinite)(p)) | |
106 { | |
107 *result = policies::raise_domain_error<RealType>( | |
108 function, | |
109 "Probability argument is %1%, but must be >= 0 and <= 1 !", p, pol); | |
110 return false; | |
111 } | |
112 return true; | |
113 } // bool check_prob | |
114 | |
115 template <class RealType, class Policy> | |
116 inline bool check_x(const char* function, const RealType& x_min, const RealType& x_max, const RealType& x, RealType* result, const Policy& pol) | |
117 { // Check x finite and x_min < x < x_max. | |
118 if (!(boost::math::isfinite)(x)) | |
119 { | |
120 *result = policies::raise_domain_error<RealType>( | |
121 function, | |
122 "x argument is %1%, but must be finite !", x, pol); | |
123 return false; | |
124 } | |
125 if ((x < x_min) || (x > x_max)) | |
126 { | |
127 // std::cout << x_min << ' ' << x << x_max << std::endl; | |
128 *result = policies::raise_domain_error<RealType>( | |
129 function, | |
130 "x argument is %1%, but must be x_min < x < x_max !", x, pol); | |
131 // For example: | |
132 // Error in function boost::math::pdf(arcsine_distribution<double> const&, double) : x argument is -1.01, but must be x_min < x < x_max ! | |
133 // TODO Perhaps show values of x_min and x_max? | |
134 return false; | |
135 } | |
136 return true; | |
137 } // bool check_x | |
138 | |
139 template <class RealType, class Policy> | |
140 inline bool check_dist(const char* function, const RealType& x_min, const RealType& x_max, RealType* result, const Policy& pol) | |
141 { // Check both x_min and x_max finite, and x_min < x_max. | |
142 return check_x_min(function, x_min, result, pol) | |
143 && check_x_max(function, x_max, result, pol) | |
144 && check_x_minmax(function, x_min, x_max, result, pol); | |
145 } // bool check_dist | |
146 | |
147 template <class RealType, class Policy> | |
148 inline bool check_dist_and_x(const char* function, const RealType& x_min, const RealType& x_max, RealType x, RealType* result, const Policy& pol) | |
149 { | |
150 return check_dist(function, x_min, x_max, result, pol) | |
151 && arcsine_detail::check_x(function, x_min, x_max, x, result, pol); | |
152 } // bool check_dist_and_x | |
153 | |
154 template <class RealType, class Policy> | |
155 inline bool check_dist_and_prob(const char* function, const RealType& x_min, const RealType& x_max, RealType p, RealType* result, const Policy& pol) | |
156 { | |
157 return check_dist(function, x_min, x_max, result, pol) | |
158 && check_prob(function, p, result, pol); | |
159 } // bool check_dist_and_prob | |
160 | |
161 } // namespace arcsine_detail | |
162 | |
163 template <class RealType = double, class Policy = policies::policy<> > | |
164 class arcsine_distribution | |
165 { | |
166 public: | |
167 typedef RealType value_type; | |
168 typedef Policy policy_type; | |
169 | |
170 arcsine_distribution(RealType x_min = 0, RealType x_max = 1) : m_x_min(x_min), m_x_max(x_max) | |
171 { // Default beta (alpha = beta = 0.5) is standard arcsine with x_min = 0, x_max = 1. | |
172 // Generalized to allow x_min and x_max to be specified. | |
173 RealType result; | |
174 arcsine_detail::check_dist( | |
175 "boost::math::arcsine_distribution<%1%>::arcsine_distribution", | |
176 m_x_min, | |
177 m_x_max, | |
178 &result, Policy()); | |
179 } // arcsine_distribution constructor. | |
180 // Accessor functions: | |
181 RealType x_min() const | |
182 { | |
183 return m_x_min; | |
184 } | |
185 RealType x_max() const | |
186 { | |
187 return m_x_max; | |
188 } | |
189 | |
190 private: | |
191 RealType m_x_min; // Two x min and x max parameters of the arcsine distribution. | |
192 RealType m_x_max; | |
193 }; // template <class RealType, class Policy> class arcsine_distribution | |
194 | |
195 // Convenient typedef to construct double version. | |
196 typedef arcsine_distribution<double> arcsine; | |
197 | |
198 | |
199 template <class RealType, class Policy> | |
200 inline const std::pair<RealType, RealType> range(const arcsine_distribution<RealType, Policy>& dist) | |
201 { // Range of permissible values for random variable x. | |
202 using boost::math::tools::max_value; | |
203 return std::pair<RealType, RealType>(static_cast<RealType>(dist.x_min()), static_cast<RealType>(dist.x_max())); | |
204 } | |
205 | |
206 template <class RealType, class Policy> | |
207 inline const std::pair<RealType, RealType> support(const arcsine_distribution<RealType, Policy>& dist) | |
208 { // Range of supported values for random variable x. | |
209 // This is range where cdf rises from 0 to 1, and outside it, the pdf is zero. | |
210 return std::pair<RealType, RealType>(static_cast<RealType>(dist.x_min()), static_cast<RealType>(dist.x_max())); | |
211 } | |
212 | |
213 template <class RealType, class Policy> | |
214 inline RealType mean(const arcsine_distribution<RealType, Policy>& dist) | |
215 { // Mean of arcsine distribution . | |
216 RealType result; | |
217 RealType x_min = dist.x_min(); | |
218 RealType x_max = dist.x_max(); | |
219 | |
220 if (false == arcsine_detail::check_dist( | |
221 "boost::math::mean(arcsine_distribution<%1%> const&, %1% )", | |
222 x_min, | |
223 x_max, | |
224 &result, Policy()) | |
225 ) | |
226 { | |
227 return result; | |
228 } | |
229 return (x_min + x_max) / 2; | |
230 } // mean | |
231 | |
232 template <class RealType, class Policy> | |
233 inline RealType variance(const arcsine_distribution<RealType, Policy>& dist) | |
234 { // Variance of standard arcsine distribution = (1-0)/8 = 0.125. | |
235 RealType result; | |
236 RealType x_min = dist.x_min(); | |
237 RealType x_max = dist.x_max(); | |
238 if (false == arcsine_detail::check_dist( | |
239 "boost::math::variance(arcsine_distribution<%1%> const&, %1% )", | |
240 x_min, | |
241 x_max, | |
242 &result, Policy()) | |
243 ) | |
244 { | |
245 return result; | |
246 } | |
247 return (x_max - x_min) * (x_max - x_min) / 8; | |
248 } // variance | |
249 | |
250 template <class RealType, class Policy> | |
251 inline RealType mode(const arcsine_distribution<RealType, Policy>& /* dist */) | |
252 { //There are always [*two] values for the mode, at ['x_min] and at ['x_max], default 0 and 1, | |
253 // so instead we raise the exception domain_error. | |
254 return policies::raise_domain_error<RealType>( | |
255 "boost::math::mode(arcsine_distribution<%1%>&)", | |
256 "The arcsine distribution has two modes at x_min and x_max: " | |
257 "so the return value is %1%.", | |
258 std::numeric_limits<RealType>::quiet_NaN(), Policy()); | |
259 } // mode | |
260 | |
261 template <class RealType, class Policy> | |
262 inline RealType median(const arcsine_distribution<RealType, Policy>& dist) | |
263 { // Median of arcsine distribution (a + b) / 2 == mean. | |
264 RealType x_min = dist.x_min(); | |
265 RealType x_max = dist.x_max(); | |
266 RealType result; | |
267 if (false == arcsine_detail::check_dist( | |
268 "boost::math::median(arcsine_distribution<%1%> const&, %1% )", | |
269 x_min, | |
270 x_max, | |
271 &result, Policy()) | |
272 ) | |
273 { | |
274 return result; | |
275 } | |
276 return (x_min + x_max) / 2; | |
277 } | |
278 | |
279 template <class RealType, class Policy> | |
280 inline RealType skewness(const arcsine_distribution<RealType, Policy>& dist) | |
281 { | |
282 RealType result; | |
283 RealType x_min = dist.x_min(); | |
284 RealType x_max = dist.x_max(); | |
285 | |
286 if (false == arcsine_detail::check_dist( | |
287 "boost::math::skewness(arcsine_distribution<%1%> const&, %1% )", | |
288 x_min, | |
289 x_max, | |
290 &result, Policy()) | |
291 ) | |
292 { | |
293 return result; | |
294 } | |
295 return 0; | |
296 } // skewness | |
297 | |
298 template <class RealType, class Policy> | |
299 inline RealType kurtosis_excess(const arcsine_distribution<RealType, Policy>& dist) | |
300 { | |
301 RealType result; | |
302 RealType x_min = dist.x_min(); | |
303 RealType x_max = dist.x_max(); | |
304 | |
305 if (false == arcsine_detail::check_dist( | |
306 "boost::math::kurtosis_excess(arcsine_distribution<%1%> const&, %1% )", | |
307 x_min, | |
308 x_max, | |
309 &result, Policy()) | |
310 ) | |
311 { | |
312 return result; | |
313 } | |
314 result = -3; | |
315 return result / 2; | |
316 } // kurtosis_excess | |
317 | |
318 template <class RealType, class Policy> | |
319 inline RealType kurtosis(const arcsine_distribution<RealType, Policy>& dist) | |
320 { | |
321 RealType result; | |
322 RealType x_min = dist.x_min(); | |
323 RealType x_max = dist.x_max(); | |
324 | |
325 if (false == arcsine_detail::check_dist( | |
326 "boost::math::kurtosis(arcsine_distribution<%1%> const&, %1% )", | |
327 x_min, | |
328 x_max, | |
329 &result, Policy()) | |
330 ) | |
331 { | |
332 return result; | |
333 } | |
334 | |
335 return 3 + kurtosis_excess(dist); | |
336 } // kurtosis | |
337 | |
338 template <class RealType, class Policy> | |
339 inline RealType pdf(const arcsine_distribution<RealType, Policy>& dist, const RealType& xx) | |
340 { // Probability Density/Mass Function arcsine. | |
341 BOOST_FPU_EXCEPTION_GUARD | |
342 BOOST_MATH_STD_USING // For ADL of std functions. | |
343 | |
344 static const char* function = "boost::math::pdf(arcsine_distribution<%1%> const&, %1%)"; | |
345 | |
346 RealType lo = dist.x_min(); | |
347 RealType hi = dist.x_max(); | |
348 RealType x = xx; | |
349 | |
350 // Argument checks: | |
351 RealType result = 0; | |
352 if (false == arcsine_detail::check_dist_and_x( | |
353 function, | |
354 lo, hi, x, | |
355 &result, Policy())) | |
356 { | |
357 return result; | |
358 } | |
359 using boost::math::constants::pi; | |
360 result = static_cast<RealType>(1) / (pi<RealType>() * sqrt((x - lo) * (hi - x))); | |
361 return result; | |
362 } // pdf | |
363 | |
364 template <class RealType, class Policy> | |
365 inline RealType cdf(const arcsine_distribution<RealType, Policy>& dist, const RealType& x) | |
366 { // Cumulative Distribution Function arcsine. | |
367 BOOST_MATH_STD_USING // For ADL of std functions. | |
368 | |
369 static const char* function = "boost::math::cdf(arcsine_distribution<%1%> const&, %1%)"; | |
370 | |
371 RealType x_min = dist.x_min(); | |
372 RealType x_max = dist.x_max(); | |
373 | |
374 // Argument checks: | |
375 RealType result = 0; | |
376 if (false == arcsine_detail::check_dist_and_x( | |
377 function, | |
378 x_min, x_max, x, | |
379 &result, Policy())) | |
380 { | |
381 return result; | |
382 } | |
383 // Special cases: | |
384 if (x == x_min) | |
385 { | |
386 return 0; | |
387 } | |
388 else if (x == x_max) | |
389 { | |
390 return 1; | |
391 } | |
392 using boost::math::constants::pi; | |
393 result = static_cast<RealType>(2) * asin(sqrt((x - x_min) / (x_max - x_min))) / pi<RealType>(); | |
394 return result; | |
395 } // arcsine cdf | |
396 | |
397 template <class RealType, class Policy> | |
398 inline RealType cdf(const complemented2_type<arcsine_distribution<RealType, Policy>, RealType>& c) | |
399 { // Complemented Cumulative Distribution Function arcsine. | |
400 BOOST_MATH_STD_USING // For ADL of std functions. | |
401 static const char* function = "boost::math::cdf(arcsine_distribution<%1%> const&, %1%)"; | |
402 | |
403 RealType x = c.param; | |
404 arcsine_distribution<RealType, Policy> const& dist = c.dist; | |
405 RealType x_min = dist.x_min(); | |
406 RealType x_max = dist.x_max(); | |
407 | |
408 // Argument checks: | |
409 RealType result = 0; | |
410 if (false == arcsine_detail::check_dist_and_x( | |
411 function, | |
412 x_min, x_max, x, | |
413 &result, Policy())) | |
414 { | |
415 return result; | |
416 } | |
417 if (x == x_min) | |
418 { | |
419 return 0; | |
420 } | |
421 else if (x == x_max) | |
422 { | |
423 return 1; | |
424 } | |
425 using boost::math::constants::pi; | |
426 // Naive version x = 1 - x; | |
427 // result = static_cast<RealType>(2) * asin(sqrt((x - x_min) / (x_max - x_min))) / pi<RealType>(); | |
428 // is less accurate, so use acos instead of asin for complement. | |
429 result = static_cast<RealType>(2) * acos(sqrt((x - x_min) / (x_max - x_min))) / pi<RealType>(); | |
430 return result; | |
431 } // arcine ccdf | |
432 | |
433 template <class RealType, class Policy> | |
434 inline RealType quantile(const arcsine_distribution<RealType, Policy>& dist, const RealType& p) | |
435 { | |
436 // Quantile or Percent Point arcsine function or | |
437 // Inverse Cumulative probability distribution function CDF. | |
438 // Return x (0 <= x <= 1), | |
439 // for a given probability p (0 <= p <= 1). | |
440 // These functions take a probability as an argument | |
441 // and return a value such that the probability that a random variable x | |
442 // will be less than or equal to that value | |
443 // is whatever probability you supplied as an argument. | |
444 BOOST_MATH_STD_USING // For ADL of std functions. | |
445 | |
446 using boost::math::constants::half_pi; | |
447 | |
448 static const char* function = "boost::math::quantile(arcsine_distribution<%1%> const&, %1%)"; | |
449 | |
450 RealType result = 0; // of argument checks: | |
451 RealType x_min = dist.x_min(); | |
452 RealType x_max = dist.x_max(); | |
453 if (false == arcsine_detail::check_dist_and_prob( | |
454 function, | |
455 x_min, x_max, p, | |
456 &result, Policy())) | |
457 { | |
458 return result; | |
459 } | |
460 // Special cases: | |
461 if (p == 0) | |
462 { | |
463 return 0; | |
464 } | |
465 if (p == 1) | |
466 { | |
467 return 1; | |
468 } | |
469 | |
470 RealType sin2hpip = sin(half_pi<RealType>() * p); | |
471 RealType sin2hpip2 = sin2hpip * sin2hpip; | |
472 result = -x_min * sin2hpip2 + x_min + x_max * sin2hpip2; | |
473 | |
474 return result; | |
475 } // quantile | |
476 | |
477 template <class RealType, class Policy> | |
478 inline RealType quantile(const complemented2_type<arcsine_distribution<RealType, Policy>, RealType>& c) | |
479 { | |
480 // Complement Quantile or Percent Point arcsine function. | |
481 // Return the number of expected x for a given | |
482 // complement of the probability q. | |
483 BOOST_MATH_STD_USING // For ADL of std functions. | |
484 | |
485 using boost::math::constants::half_pi; | |
486 static const char* function = "boost::math::quantile(arcsine_distribution<%1%> const&, %1%)"; | |
487 | |
488 // Error checks: | |
489 RealType q = c.param; | |
490 const arcsine_distribution<RealType, Policy>& dist = c.dist; | |
491 RealType result = 0; | |
492 RealType x_min = dist.x_min(); | |
493 RealType x_max = dist.x_max(); | |
494 if (false == arcsine_detail::check_dist_and_prob( | |
495 function, | |
496 x_min, | |
497 x_max, | |
498 q, | |
499 &result, Policy())) | |
500 { | |
501 return result; | |
502 } | |
503 // Special cases: | |
504 if (q == 1) | |
505 { | |
506 return 0; | |
507 } | |
508 if (q == 0) | |
509 { | |
510 return 1; | |
511 } | |
512 // Naive RealType p = 1 - q; result = sin(half_pi<RealType>() * p); loses accuracy, so use a cos alternative instead. | |
513 //result = cos(half_pi<RealType>() * q); // for arcsine(0,1) | |
514 //result = result * result; | |
515 // For generalized arcsine: | |
516 RealType cos2hpip = cos(half_pi<RealType>() * q); | |
517 RealType cos2hpip2 = cos2hpip * cos2hpip; | |
518 result = -x_min * cos2hpip2 + x_min + x_max * cos2hpip2; | |
519 | |
520 return result; | |
521 } // Quantile Complement | |
522 | |
523 } // namespace math | |
524 } // namespace boost | |
525 | |
526 // This include must be at the end, *after* the accessors | |
527 // for this distribution have been defined, in order to | |
528 // keep compilers that support two-phase lookup happy. | |
529 #include <boost/math/distributions/detail/derived_accessors.hpp> | |
530 | |
531 #if defined (BOOST_MSVC) | |
532 # pragma warning(pop) | |
533 #endif | |
534 | |
535 #endif // BOOST_MATH_DIST_ARCSINE_HPP |