]>
Commit | Line | Data |
---|---|---|
7c673cae FG |
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 | ||
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 |