math/distributions/arcsine.hpp

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