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