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> check_x_min(const char * function,const RealType & x,RealType * result,const Policy & pol)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> check_x_max(const char * function,const RealType & x,RealType * result,const Policy & pol)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> check_x_minmax(const char * function,const RealType & x_min,const RealType & x_max,RealType * result,const Policy & pol)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> check_prob(const char * function,const RealType & p,RealType * result,const Policy & pol)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> check_x(const char * function,const RealType & x_min,const RealType & x_max,const RealType & x,RealType * result,const Policy & pol)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> check_dist(const char * function,const RealType & x_min,const RealType & x_max,RealType * result,const Policy & pol)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> check_dist_and_x(const char * function,const RealType & x_min,const RealType & x_max,RealType x,RealType * result,const Policy & pol)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> check_dist_and_prob(const char * function,const RealType & x_min,const RealType & x_max,RealType p,RealType * result,const Policy & pol)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 arcsine_distribution(RealType x_min=0,RealType x_max=1)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: x_min() const181 RealType x_min() const 182 { 183 return m_x_min; 184 } x_max() const185 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> range(const arcsine_distribution<RealType,Policy> & dist)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> support(const arcsine_distribution<RealType,Policy> & dist)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> mean(const arcsine_distribution<RealType,Policy> & dist)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> variance(const arcsine_distribution<RealType,Policy> & dist)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> mode(const arcsine_distribution<RealType,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> median(const arcsine_distribution<RealType,Policy> & dist)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> skewness(const arcsine_distribution<RealType,Policy> & dist)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> kurtosis_excess(const arcsine_distribution<RealType,Policy> & dist)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> kurtosis(const arcsine_distribution<RealType,Policy> & dist)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> pdf(const arcsine_distribution<RealType,Policy> & dist,const RealType & xx)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> cdf(const arcsine_distribution<RealType,Policy> & dist,const RealType & x)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> cdf(const complemented2_type<arcsine_distribution<RealType,Policy>,RealType> & c)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> quantile(const arcsine_distribution<RealType,Policy> & dist,const RealType & p)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> quantile(const complemented2_type<arcsine_distribution<RealType,Policy>,RealType> & c)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 536