math/distributions/arcsine.hpp

0001 // boost/math/distributions/arcsine.hpp
0002
0003 // Copyright John Maddock 2014.
0004 // Copyright Paul A. Bristow 2014.
0005 // Copyright Matt Borland 2024.
0006
0007 // Use, modification and distribution are subject to the
0008 // Boost Software License, Version 1.0.
0009 // (See accompanying file LICENSE_1_0.txt
0010 // or copy at http://www.boost.org/LICENSE_1_0.txt)
0011
0012 // http://en.wikipedia.org/wiki/arcsine_distribution
0013
0014 // The arcsine Distribution is a continuous probability distribution.
0015 // http://en.wikipedia.org/wiki/Arcsine_distribution
0016 // http://www.wolframalpha.com/input/?i=ArcSinDistribution
0017
0018 // Standard arcsine distribution is a special case of beta distribution with both a & b = one half,
0019 // and 0 <= x <= 1.
0020
0021 // It is generalized to include any bounded support a <= x <= b from 0 <= x <= 1
0022 // by Wolfram and Wikipedia,
0023 // but using location and scale parameters by
0024 // Virtual Laboratories in Probability and Statistics http://www.math.uah.edu/stat/index.html
0025 // http://www.math.uah.edu/stat/special/Arcsine.html
0026 // The end-point version is simpler and more obvious, so we implement that.
0027 // TODO Perhaps provide location and scale functions?
0028
0029
0030 #ifndef BOOST_MATH_DIST_ARCSINE_HPP
0031 #define BOOST_MATH_DIST_ARCSINE_HPP
0032
0033 #include <boost/math/tools/config.hpp>
0034 #include <boost/math/tools/tuple.hpp>
0035 #include <boost/math/distributions/complement.hpp> // complements.
0036 #include <boost/math/distributions/detail/common_error_handling.hpp> // error checks.
0037 #include <boost/math/constants/constants.hpp>
0038 #include <boost/math/special_functions/fpclassify.hpp> // isnan.
0039 #include <boost/math/policies/policy.hpp>
0040 #include <boost/math/policies/error_handling.hpp>
0041
0042 #ifndef BOOST_MATH_HAS_NVRTC
0043 #include <boost/math/distributions/fwd.hpp>
0044 #include <cmath>
0045 #include <utility>
0046 #include <exception>  // For std::domain_error.
0047 #endif
0048
0049 #if defined (BOOST_MSVC)
0050 #  pragma warning(push)
0051 #  pragma warning(disable: 4702) // Unreachable code,
0052 // in domain_error_imp in error_handling.
0053 #endif
0054
0055 namespace boost
0056 {
0057   namespace math
0058   {
0059     namespace arcsine_detail
0060     {
0061       // Common error checking routines for arcsine distribution functions:
0062       // Duplicating for x_min and x_max provides specific error messages.
0063       template <class RealType, class Policy>
0064       BOOST_MATH_GPU_ENABLED inline bool check_x_min(const char* function, const RealType& x, RealType* result, const Policy& pol)
0065       {
0066         if (!(boost::math::isfinite)(x))
0067         {
0068           *result = policies::raise_domain_error<RealType>(
0069             function,
0070             "x_min argument is %1%, but must be finite !", x, pol);
0071           return false;
0072         }
0073         return true;
0074       } // bool check_x_min
0075
0076       template <class RealType, class Policy>
0077       BOOST_MATH_GPU_ENABLED inline bool check_x_max(const char* function, const RealType& x, RealType* result, const Policy& pol)
0078       {
0079         if (!(boost::math::isfinite)(x))
0080         {
0081           *result = policies::raise_domain_error<RealType>(
0082             function,
0083             "x_max argument is %1%, but must be finite !", x, pol);
0084           return false;
0085         }
0086         return true;
0087       } // bool check_x_max
0088
0089
0090       template <class RealType, class Policy>
0091       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)
0092       { // Check x_min < x_max
0093         if (x_min >= x_max)
0094         {
0095           constexpr auto msg = "x_max argument is %1%, but must be > x_min";
0096           *result = policies::raise_domain_error<RealType>(
0097             function,
0098             msg, x_max, pol);
0099             // "x_max argument is %1%, but must be > x_min !", x_max, pol);
0100             // "x_max argument is %1%, but must be > x_min %2!", x_max, x_min, pol); would be better.
0101             // But would require replication of all helpers functions in /policies/error_handling.hpp for two values,
0102             // as well as two value versions of raise_error, raise_domain_error and do_format
0103           return false;
0104         }
0105         return true;
0106       } // bool check_x_minmax
0107
0108       template <class RealType, class Policy>
0109       BOOST_MATH_GPU_ENABLED inline bool check_prob(const char* function, const RealType& p, RealType* result, const Policy& pol)
0110       {
0111         if ((p < 0) || (p > 1) || !(boost::math::isfinite)(p))
0112         {
0113           *result = policies::raise_domain_error<RealType>(
0114             function,
0115             "Probability argument is %1%, but must be >= 0 and <= 1 !", p, pol);
0116           return false;
0117         }
0118         return true;
0119       } // bool check_prob
0120
0121       template <class RealType, class Policy>
0122       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)
0123       { // Check x finite and x_min < x < x_max.
0124         if (!(boost::math::isfinite)(x))
0125         {
0126           *result = policies::raise_domain_error<RealType>(
0127             function,
0128             "x argument is %1%, but must be finite !", x, pol);
0129           return false;
0130         }
0131         if ((x < x_min) || (x > x_max))
0132         {
0133           // std::cout << x_min << ' ' << x << x_max << std::endl;
0134           *result = policies::raise_domain_error<RealType>(
0135             function,
0136             "x argument is %1%, but must be x_min < x < x_max !", x, pol);
0137           // For example:
0138           // Error in function boost::math::pdf(arcsine_distribution<double> const&, double) : x argument is -1.01, but must be x_min < x < x_max !
0139           // TODO Perhaps show values of x_min and x_max?
0140           return false;
0141         }
0142         return true;
0143       } // bool check_x
0144
0145       template <class RealType, class Policy>
0146       BOOST_MATH_GPU_ENABLED inline bool check_dist(const char* function, const RealType& x_min, const RealType& x_max, RealType* result, const Policy& pol)
0147       { // Check both x_min and x_max finite, and x_min  < x_max.
0148         return check_x_min(function, x_min, result, pol)
0149             && check_x_max(function, x_max, result, pol)
0150             && check_x_minmax(function, x_min, x_max, result, pol);
0151       } // bool check_dist
0152
0153       template <class RealType, class Policy>
0154       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)
0155       {
0156         return check_dist(function, x_min, x_max, result, pol)
0157           && arcsine_detail::check_x(function, x_min, x_max, x, result, pol);
0158       } // bool check_dist_and_x
0159
0160       template <class RealType, class Policy>
0161       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)
0162       {
0163         return check_dist(function, x_min, x_max, result, pol)
0164           && check_prob(function, p, result, pol);
0165       } // bool check_dist_and_prob
0166
0167     } // namespace arcsine_detail
0168
0169     template <class RealType = double, class Policy = policies::policy<> >
0170     class arcsine_distribution
0171     {
0172     public:
0173       typedef RealType value_type;
0174       typedef Policy policy_type;
0175
0176       BOOST_MATH_GPU_ENABLED arcsine_distribution(RealType x_min = 0, RealType x_max = 1) : m_x_min(x_min), m_x_max(x_max)
0177       { // Default beta (alpha = beta = 0.5) is standard arcsine with x_min = 0, x_max = 1.
0178         // Generalized to allow x_min and x_max to be specified.
0179         RealType result;
0180         arcsine_detail::check_dist(
0181           "boost::math::arcsine_distribution<%1%>::arcsine_distribution",
0182           m_x_min,
0183           m_x_max,
0184           &result, Policy());
0185       } // arcsine_distribution constructor.
0186       // Accessor functions:
0187       BOOST_MATH_GPU_ENABLED RealType x_min() const
0188       {
0189         return m_x_min;
0190       }
0191       BOOST_MATH_GPU_ENABLED RealType x_max() const
0192       {
0193         return m_x_max;
0194       }
0195
0196     private:
0197       RealType m_x_min; // Two x min and x max parameters of the arcsine distribution.
0198       RealType m_x_max;
0199     }; // template <class RealType, class Policy> class arcsine_distribution
0200
0201     // Convenient typedef to construct double version.
0202     typedef arcsine_distribution<double> arcsine;
0203
0204     #ifdef __cpp_deduction_guides
0205     template <class RealType>
0206     arcsine_distribution(RealType)->arcsine_distribution<typename boost::math::tools::promote_args<RealType>::type>;
0207     template <class RealType>
0208     arcsine_distribution(RealType, RealType)->arcsine_distribution<typename boost::math::tools::promote_args<RealType>::type>;
0209     #endif
0210
0211     template <class RealType, class Policy>
0212     BOOST_MATH_GPU_ENABLED inline const boost::math::pair<RealType, RealType> range(const arcsine_distribution<RealType, Policy>&  dist)
0213     { // Range of permissible values for random variable x.
0214       using boost::math::tools::max_value;
0215       return boost::math::pair<RealType, RealType>(static_cast<RealType>(dist.x_min()), static_cast<RealType>(dist.x_max()));
0216     }
0217
0218     template <class RealType, class Policy>
0219     BOOST_MATH_GPU_ENABLED inline const boost::math::pair<RealType, RealType> support(const arcsine_distribution<RealType, Policy>&  dist)
0220     { // Range of supported values for random variable x.
0221       // This is range where cdf rises from 0 to 1, and outside it, the pdf is zero.
0222       return boost::math::pair<RealType, RealType>(static_cast<RealType>(dist.x_min()), static_cast<RealType>(dist.x_max()));
0223     }
0224
0225     template <class RealType, class Policy>
0226     BOOST_MATH_GPU_ENABLED inline RealType mean(const arcsine_distribution<RealType, Policy>& dist)
0227     { // Mean of arcsine distribution .
0228       RealType result;
0229       RealType x_min = dist.x_min();
0230       RealType x_max = dist.x_max();
0231
0232       if (false == arcsine_detail::check_dist(
0233         "boost::math::mean(arcsine_distribution<%1%> const&, %1% )",
0234         x_min,
0235         x_max,
0236         &result, Policy())
0237         )
0238       {
0239         return result;
0240       }
0241       return  (x_min + x_max) / 2;
0242     } // mean
0243
0244     template <class RealType, class Policy>
0245     BOOST_MATH_GPU_ENABLED inline RealType variance(const arcsine_distribution<RealType, Policy>& dist)
0246     { // Variance of standard arcsine distribution = (1-0)/8 = 0.125.
0247       RealType result;
0248       RealType x_min = dist.x_min();
0249       RealType x_max = dist.x_max();
0250       if (false == arcsine_detail::check_dist(
0251         "boost::math::variance(arcsine_distribution<%1%> const&, %1% )",
0252         x_min,
0253         x_max,
0254         &result, Policy())
0255         )
0256       {
0257         return result;
0258       }
0259       return  (x_max - x_min) * (x_max - x_min) / 8;
0260     } // variance
0261
0262     template <class RealType, class Policy>
0263     BOOST_MATH_GPU_ENABLED inline RealType mode(const arcsine_distribution<RealType, Policy>& /* dist */)
0264     { //There are always [*two] values for the mode, at ['x_min] and at ['x_max], default 0 and 1,
0265       // so instead we raise the exception domain_error.
0266       return policies::raise_domain_error<RealType>(
0267         "boost::math::mode(arcsine_distribution<%1%>&)",
0268         "The arcsine distribution has two modes at x_min and x_max: "
0269         "so the return value is %1%.",
0270         std::numeric_limits<RealType>::quiet_NaN(), Policy());
0271     } // mode
0272
0273     template <class RealType, class Policy>
0274     BOOST_MATH_GPU_ENABLED inline RealType median(const arcsine_distribution<RealType, Policy>& dist)
0275     { // Median of arcsine distribution (a + b) / 2 == mean.
0276       RealType x_min = dist.x_min();
0277       RealType x_max = dist.x_max();
0278       RealType result;
0279       if (false == arcsine_detail::check_dist(
0280         "boost::math::median(arcsine_distribution<%1%> const&, %1% )",
0281         x_min,
0282         x_max,
0283         &result, Policy())
0284         )
0285       {
0286         return result;
0287       }
0288       return  (x_min + x_max) / 2;
0289     }
0290
0291     template <class RealType, class Policy>
0292     BOOST_MATH_GPU_ENABLED inline RealType skewness(const arcsine_distribution<RealType, Policy>& dist)
0293     {
0294       RealType result;
0295       RealType x_min = dist.x_min();
0296       RealType x_max = dist.x_max();
0297
0298       if (false == arcsine_detail::check_dist(
0299         "boost::math::skewness(arcsine_distribution<%1%> const&, %1% )",
0300         x_min,
0301         x_max,
0302         &result, Policy())
0303         )
0304       {
0305         return result;
0306       }
0307       return 0;
0308     } // skewness
0309
0310     template <class RealType, class Policy>
0311     BOOST_MATH_GPU_ENABLED inline RealType kurtosis_excess(const arcsine_distribution<RealType, Policy>& dist)
0312     {
0313       RealType result;
0314       RealType x_min = dist.x_min();
0315       RealType x_max = dist.x_max();
0316
0317       if (false == arcsine_detail::check_dist(
0318         "boost::math::kurtosis_excess(arcsine_distribution<%1%> const&, %1% )",
0319         x_min,
0320         x_max,
0321         &result, Policy())
0322         )
0323       {
0324         return result;
0325       }
0326       result = -3;
0327       return  result / 2;
0328     } // kurtosis_excess
0329
0330     template <class RealType, class Policy>
0331     BOOST_MATH_GPU_ENABLED inline RealType kurtosis(const arcsine_distribution<RealType, Policy>& dist)
0332     {
0333       RealType result;
0334       RealType x_min = dist.x_min();
0335       RealType x_max = dist.x_max();
0336
0337       if (false == arcsine_detail::check_dist(
0338         "boost::math::kurtosis(arcsine_distribution<%1%> const&, %1% )",
0339         x_min,
0340         x_max,
0341         &result, Policy())
0342         )
0343       {
0344         return result;
0345       }
0346
0347       return 3 + kurtosis_excess(dist);
0348     } // kurtosis
0349
0350     template <class RealType, class Policy>
0351     BOOST_MATH_GPU_ENABLED inline RealType pdf(const arcsine_distribution<RealType, Policy>& dist, const RealType& xx)
0352     { // Probability Density/Mass Function arcsine.
0353       BOOST_FPU_EXCEPTION_GUARD
0354       BOOST_MATH_STD_USING // For ADL of std functions.
0355
0356       constexpr auto function = "boost::math::pdf(arcsine_distribution<%1%> const&, %1%)";
0357
0358       RealType lo = dist.x_min();
0359       RealType hi = dist.x_max();
0360       RealType x = xx;
0361
0362       // Argument checks:
0363       RealType result = 0;
0364       if (false == arcsine_detail::check_dist_and_x(
0365         function,
0366         lo, hi, x,
0367         &result, Policy()))
0368       {
0369         return result;
0370       }
0371       using boost::math::constants::pi;
0372       result = static_cast<RealType>(1) / (pi<RealType>() * sqrt((x - lo) * (hi - x)));
0373       return result;
0374     } // pdf
0375
0376     template <class RealType, class Policy>
0377     BOOST_MATH_GPU_ENABLED inline RealType cdf(const arcsine_distribution<RealType, Policy>& dist, const RealType& x)
0378     { // Cumulative Distribution Function arcsine.
0379       BOOST_MATH_STD_USING // For ADL of std functions.
0380
0381       constexpr auto function = "boost::math::cdf(arcsine_distribution<%1%> const&, %1%)";
0382
0383       RealType x_min = dist.x_min();
0384       RealType x_max = dist.x_max();
0385
0386       // Argument checks:
0387       RealType result = 0;
0388       if (false == arcsine_detail::check_dist_and_x(
0389         function,
0390         x_min, x_max, x,
0391         &result, Policy()))
0392       {
0393         return result;
0394       }
0395       // Special cases:
0396       if (x == x_min)
0397       {
0398         return 0;
0399       }
0400       else if (x == x_max)
0401       {
0402         return 1;
0403       }
0404       using boost::math::constants::pi;
0405       result = static_cast<RealType>(2) * asin(sqrt((x - x_min) / (x_max - x_min))) / pi<RealType>();
0406       return result;
0407     } // arcsine cdf
0408
0409     template <class RealType, class Policy>
0410     BOOST_MATH_GPU_ENABLED inline RealType cdf(const complemented2_type<arcsine_distribution<RealType, Policy>, RealType>& c)
0411     { // Complemented Cumulative Distribution Function arcsine.
0412       BOOST_MATH_STD_USING // For ADL of std functions.
0413       constexpr auto function = "boost::math::cdf(arcsine_distribution<%1%> const&, %1%)";
0414
0415       RealType x = c.param;
0416       arcsine_distribution<RealType, Policy> const& dist = c.dist;
0417       RealType x_min = dist.x_min();
0418       RealType x_max = dist.x_max();
0419
0420       // Argument checks:
0421       RealType result = 0;
0422       if (false == arcsine_detail::check_dist_and_x(
0423         function,
0424         x_min, x_max, x,
0425         &result, Policy()))
0426       {
0427         return result;
0428       }
0429       if (x == x_min)
0430       {
0431         return 0;
0432       }
0433       else if (x == x_max)
0434       {
0435         return 1;
0436       }
0437       using boost::math::constants::pi;
0438       // Naive version x = 1 - x;
0439       // result = static_cast<RealType>(2) * asin(sqrt((x - x_min) / (x_max - x_min))) / pi<RealType>();
0440       // is less accurate, so use acos instead of asin for complement.
0441       result = static_cast<RealType>(2) * acos(sqrt((x - x_min) / (x_max - x_min))) / pi<RealType>();
0442       return result;
0443     } // arcsine ccdf
0444
0445     template <class RealType, class Policy>
0446     BOOST_MATH_GPU_ENABLED inline RealType quantile(const arcsine_distribution<RealType, Policy>& dist, const RealType& p)
0447     {
0448       // Quantile or Percent Point arcsine function or
0449       // Inverse Cumulative probability distribution function CDF.
0450       // Return x (0 <= x <= 1),
0451       // for a given probability p (0 <= p <= 1).
0452       // These functions take a probability as an argument
0453       // and return a value such that the probability that a random variable x
0454       // will be less than or equal to that value
0455       // is whatever probability you supplied as an argument.
0456       BOOST_MATH_STD_USING // For ADL of std functions.
0457
0458       using boost::math::constants::half_pi;
0459
0460       constexpr auto function = "boost::math::quantile(arcsine_distribution<%1%> const&, %1%)";
0461
0462       RealType result = 0; // of argument checks:
0463       RealType x_min = dist.x_min();
0464       RealType x_max = dist.x_max();
0465       if (false == arcsine_detail::check_dist_and_prob(
0466         function,
0467         x_min, x_max, p,
0468         &result, Policy()))
0469       {
0470         return result;
0471       }
0472       // Special cases:
0473       if (p == 0)
0474       {
0475         return 0;
0476       }
0477       if (p == 1)
0478       {
0479         return 1;
0480       }
0481
0482       RealType sin2hpip = sin(half_pi<RealType>() * p);
0483       RealType sin2hpip2 = sin2hpip * sin2hpip;
0484       result = -x_min * sin2hpip2 + x_min + x_max * sin2hpip2;
0485
0486       return result;
0487     } // quantile
0488
0489     template <class RealType, class Policy>
0490     BOOST_MATH_GPU_ENABLED inline RealType quantile(const complemented2_type<arcsine_distribution<RealType, Policy>, RealType>& c)
0491     {
0492       // Complement Quantile or Percent Point arcsine function.
0493       // Return the number of expected x for a given
0494       // complement of the probability q.
0495       BOOST_MATH_STD_USING // For ADL of std functions.
0496
0497       using boost::math::constants::half_pi;
0498       constexpr auto function = "boost::math::quantile(arcsine_distribution<%1%> const&, %1%)";
0499
0500       // Error checks:
0501       RealType q = c.param;
0502       const arcsine_distribution<RealType, Policy>& dist = c.dist;
0503       RealType result = 0;
0504       RealType x_min = dist.x_min();
0505       RealType x_max = dist.x_max();
0506       if (false == arcsine_detail::check_dist_and_prob(
0507         function,
0508         x_min,
0509         x_max,
0510         q,
0511         &result, Policy()))
0512       {
0513         return result;
0514       }
0515       // Special cases:
0516       if (q == 1)
0517       {
0518         return 0;
0519       }
0520       if (q == 0)
0521       {
0522         return 1;
0523       }
0524       // Naive RealType p = 1 - q; result = sin(half_pi<RealType>() * p); loses accuracy, so use a cos alternative instead.
0525       //result = cos(half_pi<RealType>() * q); // for arcsine(0,1)
0526       //result = result * result;
0527       // For generalized arcsine:
0528       RealType cos2hpip = cos(half_pi<RealType>() * q);
0529       RealType cos2hpip2 = cos2hpip * cos2hpip;
0530       result = -x_min * cos2hpip2 + x_min + x_max * cos2hpip2;
0531
0532       return result;
0533     } // Quantile Complement
0534
0535   } // namespace math
0536 } // namespace boost
0537
0538 // This include must be at the end, *after* the accessors
0539 // for this distribution have been defined, in order to
0540 // keep compilers that support two-phase lookup happy.
0541 #include <boost/math/distributions/detail/derived_accessors.hpp>
0542
0543 #if defined (BOOST_MSVC)
0544 # pragma warning(pop)
0545 #endif
0546
0547 #endif // BOOST_MATH_DIST_ARCSINE_HPP