Chris@102: // boost/math/distributions/arcsine.hpp Chris@102: Chris@102: // Copyright John Maddock 2014. Chris@102: // Copyright Paul A. Bristow 2014. Chris@102: Chris@102: // Use, modification and distribution are subject to the Chris@102: // Boost Software License, Version 1.0. Chris@102: // (See accompanying file LICENSE_1_0.txt Chris@102: // or copy at http://www.boost.org/LICENSE_1_0.txt) Chris@102: Chris@102: // http://en.wikipedia.org/wiki/arcsine_distribution Chris@102: Chris@102: // The arcsine Distribution is a continuous probability distribution. Chris@102: // http://en.wikipedia.org/wiki/Arcsine_distribution Chris@102: // http://www.wolframalpha.com/input/?i=ArcSinDistribution Chris@102: Chris@102: // Standard arcsine distribution is a special case of beta distribution with both a & b = one half, Chris@102: // and 0 <= x <= 1. Chris@102: Chris@102: // It is generalized to include any bounded support a <= x <= b from 0 <= x <= 1 Chris@102: // by Wolfram and Wikipedia, Chris@102: // but using location and scale parameters by Chris@102: // Virtual Laboratories in Probability and Statistics http://www.math.uah.edu/stat/index.html Chris@102: // http://www.math.uah.edu/stat/special/Arcsine.html Chris@102: // The end-point version is simpler and more obvious, so we implement that. Chris@102: // TODO Perhaps provide location and scale functions? Chris@102: Chris@102: Chris@102: #ifndef BOOST_MATH_DIST_ARCSINE_HPP Chris@102: #define BOOST_MATH_DIST_ARCSINE_HPP Chris@102: Chris@102: #include Chris@102: #include // complements. Chris@102: #include // error checks. Chris@102: #include Chris@102: Chris@102: #include // isnan. Chris@102: Chris@102: #if defined (BOOST_MSVC) Chris@102: # pragma warning(push) Chris@102: # pragma warning(disable: 4702) // Unreachable code, Chris@102: // in domain_error_imp in error_handling. Chris@102: #endif Chris@102: Chris@102: #include Chris@102: #include // For std::domain_error. Chris@102: Chris@102: namespace boost Chris@102: { Chris@102: namespace math Chris@102: { Chris@102: namespace arcsine_detail Chris@102: { Chris@102: // Common error checking routines for arcsine distribution functions: Chris@102: // Duplicating for x_min and x_max provides specific error messages. Chris@102: template Chris@102: inline bool check_x_min(const char* function, const RealType& x, RealType* result, const Policy& pol) Chris@102: { Chris@102: if (!(boost::math::isfinite)(x)) Chris@102: { Chris@102: *result = policies::raise_domain_error( Chris@102: function, Chris@102: "x_min argument is %1%, but must be finite !", x, pol); Chris@102: return false; Chris@102: } Chris@102: return true; Chris@102: } // bool check_x_min Chris@102: Chris@102: template Chris@102: inline bool check_x_max(const char* function, const RealType& x, RealType* result, const Policy& pol) Chris@102: { Chris@102: if (!(boost::math::isfinite)(x)) Chris@102: { Chris@102: *result = policies::raise_domain_error( Chris@102: function, Chris@102: "x_max argument is %1%, but must be finite !", x, pol); Chris@102: return false; Chris@102: } Chris@102: return true; Chris@102: } // bool check_x_max Chris@102: Chris@102: Chris@102: template Chris@102: inline bool check_x_minmax(const char* function, const RealType& x_min, const RealType& x_max, RealType* result, const Policy& pol) Chris@102: { // Check x_min < x_max Chris@102: if (x_min >= x_max) Chris@102: { Chris@102: std::string msg = "x_max argument is %1%, but must be > x_min = " + lexical_cast(x_min) + "!"; Chris@102: *result = policies::raise_domain_error( Chris@102: function, Chris@102: msg.c_str(), x_max, pol); Chris@102: // "x_max argument is %1%, but must be > x_min !", x_max, pol); Chris@102: // "x_max argument is %1%, but must be > x_min %2!", x_max, x_min, pol); would be better. Chris@102: // But would require replication of all helpers functions in /policies/error_handling.hpp for two values, Chris@102: // as well as two value versions of raise_error, raise_domain_error and do_format ... Chris@102: // so use slightly hacky lexical_cast to string instead. Chris@102: return false; Chris@102: } Chris@102: return true; Chris@102: } // bool check_x_minmax Chris@102: Chris@102: template Chris@102: inline bool check_prob(const char* function, const RealType& p, RealType* result, const Policy& pol) Chris@102: { Chris@102: if ((p < 0) || (p > 1) || !(boost::math::isfinite)(p)) Chris@102: { Chris@102: *result = policies::raise_domain_error( Chris@102: function, Chris@102: "Probability argument is %1%, but must be >= 0 and <= 1 !", p, pol); Chris@102: return false; Chris@102: } Chris@102: return true; Chris@102: } // bool check_prob Chris@102: Chris@102: template Chris@102: inline bool check_x(const char* function, const RealType& x_min, const RealType& x_max, const RealType& x, RealType* result, const Policy& pol) Chris@102: { // Check x finite and x_min < x < x_max. Chris@102: if (!(boost::math::isfinite)(x)) Chris@102: { Chris@102: *result = policies::raise_domain_error( Chris@102: function, Chris@102: "x argument is %1%, but must be finite !", x, pol); Chris@102: return false; Chris@102: } Chris@102: if ((x < x_min) || (x > x_max)) Chris@102: { Chris@102: // std::cout << x_min << ' ' << x << x_max << std::endl; Chris@102: *result = policies::raise_domain_error( Chris@102: function, Chris@102: "x argument is %1%, but must be x_min < x < x_max !", x, pol); Chris@102: // For example: Chris@102: // Error in function boost::math::pdf(arcsine_distribution const&, double) : x argument is -1.01, but must be x_min < x < x_max ! Chris@102: // TODO Perhaps show values of x_min and x_max? Chris@102: return false; Chris@102: } Chris@102: return true; Chris@102: } // bool check_x Chris@102: Chris@102: template Chris@102: inline bool check_dist(const char* function, const RealType& x_min, const RealType& x_max, RealType* result, const Policy& pol) Chris@102: { // Check both x_min and x_max finite, and x_min < x_max. Chris@102: return check_x_min(function, x_min, result, pol) Chris@102: && check_x_max(function, x_max, result, pol) Chris@102: && check_x_minmax(function, x_min, x_max, result, pol); Chris@102: } // bool check_dist Chris@102: Chris@102: template Chris@102: inline bool check_dist_and_x(const char* function, const RealType& x_min, const RealType& x_max, RealType x, RealType* result, const Policy& pol) Chris@102: { Chris@102: return check_dist(function, x_min, x_max, result, pol) Chris@102: && arcsine_detail::check_x(function, x_min, x_max, x, result, pol); Chris@102: } // bool check_dist_and_x Chris@102: Chris@102: template Chris@102: inline bool check_dist_and_prob(const char* function, const RealType& x_min, const RealType& x_max, RealType p, RealType* result, const Policy& pol) Chris@102: { Chris@102: return check_dist(function, x_min, x_max, result, pol) Chris@102: && check_prob(function, p, result, pol); Chris@102: } // bool check_dist_and_prob Chris@102: Chris@102: } // namespace arcsine_detail Chris@102: Chris@102: template > Chris@102: class arcsine_distribution Chris@102: { Chris@102: public: Chris@102: typedef RealType value_type; Chris@102: typedef Policy policy_type; Chris@102: Chris@102: arcsine_distribution(RealType x_min = 0, RealType x_max = 1) : m_x_min(x_min), m_x_max(x_max) Chris@102: { // Default beta (alpha = beta = 0.5) is standard arcsine with x_min = 0, x_max = 1. Chris@102: // Generalized to allow x_min and x_max to be specified. Chris@102: RealType result; Chris@102: arcsine_detail::check_dist( Chris@102: "boost::math::arcsine_distribution<%1%>::arcsine_distribution", Chris@102: m_x_min, Chris@102: m_x_max, Chris@102: &result, Policy()); Chris@102: } // arcsine_distribution constructor. Chris@102: // Accessor functions: Chris@102: RealType x_min() const Chris@102: { Chris@102: return m_x_min; Chris@102: } Chris@102: RealType x_max() const Chris@102: { Chris@102: return m_x_max; Chris@102: } Chris@102: Chris@102: private: Chris@102: RealType m_x_min; // Two x min and x max parameters of the arcsine distribution. Chris@102: RealType m_x_max; Chris@102: }; // template class arcsine_distribution Chris@102: Chris@102: // Convenient typedef to construct double version. Chris@102: typedef arcsine_distribution arcsine; Chris@102: Chris@102: Chris@102: template Chris@102: inline const std::pair range(const arcsine_distribution& dist) Chris@102: { // Range of permissible values for random variable x. Chris@102: using boost::math::tools::max_value; Chris@102: return std::pair(static_cast(dist.x_min()), static_cast(dist.x_max())); Chris@102: } Chris@102: Chris@102: template Chris@102: inline const std::pair support(const arcsine_distribution& dist) Chris@102: { // Range of supported values for random variable x. Chris@102: // This is range where cdf rises from 0 to 1, and outside it, the pdf is zero. Chris@102: return std::pair(static_cast(dist.x_min()), static_cast(dist.x_max())); Chris@102: } Chris@102: Chris@102: template Chris@102: inline RealType mean(const arcsine_distribution& dist) Chris@102: { // Mean of arcsine distribution . Chris@102: RealType result; Chris@102: RealType x_min = dist.x_min(); Chris@102: RealType x_max = dist.x_max(); Chris@102: Chris@102: if (false == arcsine_detail::check_dist( Chris@102: "boost::math::mean(arcsine_distribution<%1%> const&, %1% )", Chris@102: x_min, Chris@102: x_max, Chris@102: &result, Policy()) Chris@102: ) Chris@102: { Chris@102: return result; Chris@102: } Chris@102: return (x_min + x_max) / 2; Chris@102: } // mean Chris@102: Chris@102: template Chris@102: inline RealType variance(const arcsine_distribution& dist) Chris@102: { // Variance of standard arcsine distribution = (1-0)/8 = 0.125. Chris@102: RealType result; Chris@102: RealType x_min = dist.x_min(); Chris@102: RealType x_max = dist.x_max(); Chris@102: if (false == arcsine_detail::check_dist( Chris@102: "boost::math::variance(arcsine_distribution<%1%> const&, %1% )", Chris@102: x_min, Chris@102: x_max, Chris@102: &result, Policy()) Chris@102: ) Chris@102: { Chris@102: return result; Chris@102: } Chris@102: return (x_max - x_min) * (x_max - x_min) / 8; Chris@102: } // variance Chris@102: Chris@102: template Chris@102: inline RealType mode(const arcsine_distribution& /* dist */) Chris@102: { //There are always [*two] values for the mode, at ['x_min] and at ['x_max], default 0 and 1, Chris@102: // so instead we raise the exception domain_error. Chris@102: return policies::raise_domain_error( Chris@102: "boost::math::mode(arcsine_distribution<%1%>&)", Chris@102: "The arcsine distribution has two modes at x_min and x_max: " Chris@102: "so the return value is %1%.", Chris@102: std::numeric_limits::quiet_NaN(), Policy()); Chris@102: } // mode Chris@102: Chris@102: template Chris@102: inline RealType median(const arcsine_distribution& dist) Chris@102: { // Median of arcsine distribution (a + b) / 2 == mean. Chris@102: RealType x_min = dist.x_min(); Chris@102: RealType x_max = dist.x_max(); Chris@102: RealType result; Chris@102: if (false == arcsine_detail::check_dist( Chris@102: "boost::math::median(arcsine_distribution<%1%> const&, %1% )", Chris@102: x_min, Chris@102: x_max, Chris@102: &result, Policy()) Chris@102: ) Chris@102: { Chris@102: return result; Chris@102: } Chris@102: return (x_min + x_max) / 2; Chris@102: } Chris@102: Chris@102: template Chris@102: inline RealType skewness(const arcsine_distribution& dist) Chris@102: { Chris@102: RealType result; Chris@102: RealType x_min = dist.x_min(); Chris@102: RealType x_max = dist.x_max(); Chris@102: Chris@102: if (false == arcsine_detail::check_dist( Chris@102: "boost::math::skewness(arcsine_distribution<%1%> const&, %1% )", Chris@102: x_min, Chris@102: x_max, Chris@102: &result, Policy()) Chris@102: ) Chris@102: { Chris@102: return result; Chris@102: } Chris@102: return 0; Chris@102: } // skewness Chris@102: Chris@102: template Chris@102: inline RealType kurtosis_excess(const arcsine_distribution& dist) Chris@102: { Chris@102: RealType result; Chris@102: RealType x_min = dist.x_min(); Chris@102: RealType x_max = dist.x_max(); Chris@102: Chris@102: if (false == arcsine_detail::check_dist( Chris@102: "boost::math::kurtosis_excess(arcsine_distribution<%1%> const&, %1% )", Chris@102: x_min, Chris@102: x_max, Chris@102: &result, Policy()) Chris@102: ) Chris@102: { Chris@102: return result; Chris@102: } Chris@102: result = -3; Chris@102: return result / 2; Chris@102: } // kurtosis_excess Chris@102: Chris@102: template Chris@102: inline RealType kurtosis(const arcsine_distribution& dist) Chris@102: { Chris@102: RealType result; Chris@102: RealType x_min = dist.x_min(); Chris@102: RealType x_max = dist.x_max(); Chris@102: Chris@102: if (false == arcsine_detail::check_dist( Chris@102: "boost::math::kurtosis(arcsine_distribution<%1%> const&, %1% )", Chris@102: x_min, Chris@102: x_max, Chris@102: &result, Policy()) Chris@102: ) Chris@102: { Chris@102: return result; Chris@102: } Chris@102: Chris@102: return 3 + kurtosis_excess(dist); Chris@102: } // kurtosis Chris@102: Chris@102: template Chris@102: inline RealType pdf(const arcsine_distribution& dist, const RealType& xx) Chris@102: { // Probability Density/Mass Function arcsine. Chris@102: BOOST_FPU_EXCEPTION_GUARD Chris@102: BOOST_MATH_STD_USING // For ADL of std functions. Chris@102: Chris@102: static const char* function = "boost::math::pdf(arcsine_distribution<%1%> const&, %1%)"; Chris@102: Chris@102: RealType lo = dist.x_min(); Chris@102: RealType hi = dist.x_max(); Chris@102: RealType x = xx; Chris@102: Chris@102: // Argument checks: Chris@102: RealType result = 0; Chris@102: if (false == arcsine_detail::check_dist_and_x( Chris@102: function, Chris@102: lo, hi, x, Chris@102: &result, Policy())) Chris@102: { Chris@102: return result; Chris@102: } Chris@102: using boost::math::constants::pi; Chris@102: result = static_cast(1) / (pi() * sqrt((x - lo) * (hi - x))); Chris@102: return result; Chris@102: } // pdf Chris@102: Chris@102: template Chris@102: inline RealType cdf(const arcsine_distribution& dist, const RealType& x) Chris@102: { // Cumulative Distribution Function arcsine. Chris@102: BOOST_MATH_STD_USING // For ADL of std functions. Chris@102: Chris@102: static const char* function = "boost::math::cdf(arcsine_distribution<%1%> const&, %1%)"; Chris@102: Chris@102: RealType x_min = dist.x_min(); Chris@102: RealType x_max = dist.x_max(); Chris@102: Chris@102: // Argument checks: Chris@102: RealType result = 0; Chris@102: if (false == arcsine_detail::check_dist_and_x( Chris@102: function, Chris@102: x_min, x_max, x, Chris@102: &result, Policy())) Chris@102: { Chris@102: return result; Chris@102: } Chris@102: // Special cases: Chris@102: if (x == x_min) Chris@102: { Chris@102: return 0; Chris@102: } Chris@102: else if (x == x_max) Chris@102: { Chris@102: return 1; Chris@102: } Chris@102: using boost::math::constants::pi; Chris@102: result = static_cast(2) * asin(sqrt((x - x_min) / (x_max - x_min))) / pi(); Chris@102: return result; Chris@102: } // arcsine cdf Chris@102: Chris@102: template Chris@102: inline RealType cdf(const complemented2_type, RealType>& c) Chris@102: { // Complemented Cumulative Distribution Function arcsine. Chris@102: BOOST_MATH_STD_USING // For ADL of std functions. Chris@102: static const char* function = "boost::math::cdf(arcsine_distribution<%1%> const&, %1%)"; Chris@102: Chris@102: RealType x = c.param; Chris@102: arcsine_distribution const& dist = c.dist; Chris@102: RealType x_min = dist.x_min(); Chris@102: RealType x_max = dist.x_max(); Chris@102: Chris@102: // Argument checks: Chris@102: RealType result = 0; Chris@102: if (false == arcsine_detail::check_dist_and_x( Chris@102: function, Chris@102: x_min, x_max, x, Chris@102: &result, Policy())) Chris@102: { Chris@102: return result; Chris@102: } Chris@102: if (x == x_min) Chris@102: { Chris@102: return 0; Chris@102: } Chris@102: else if (x == x_max) Chris@102: { Chris@102: return 1; Chris@102: } Chris@102: using boost::math::constants::pi; Chris@102: // Naive version x = 1 - x; Chris@102: // result = static_cast(2) * asin(sqrt((x - x_min) / (x_max - x_min))) / pi(); Chris@102: // is less accurate, so use acos instead of asin for complement. Chris@102: result = static_cast(2) * acos(sqrt((x - x_min) / (x_max - x_min))) / pi(); Chris@102: return result; Chris@102: } // arcine ccdf Chris@102: Chris@102: template Chris@102: inline RealType quantile(const arcsine_distribution& dist, const RealType& p) Chris@102: { Chris@102: // Quantile or Percent Point arcsine function or Chris@102: // Inverse Cumulative probability distribution function CDF. Chris@102: // Return x (0 <= x <= 1), Chris@102: // for a given probability p (0 <= p <= 1). Chris@102: // These functions take a probability as an argument Chris@102: // and return a value such that the probability that a random variable x Chris@102: // will be less than or equal to that value Chris@102: // is whatever probability you supplied as an argument. Chris@102: BOOST_MATH_STD_USING // For ADL of std functions. Chris@102: Chris@102: using boost::math::constants::half_pi; Chris@102: Chris@102: static const char* function = "boost::math::quantile(arcsine_distribution<%1%> const&, %1%)"; Chris@102: Chris@102: RealType result = 0; // of argument checks: Chris@102: RealType x_min = dist.x_min(); Chris@102: RealType x_max = dist.x_max(); Chris@102: if (false == arcsine_detail::check_dist_and_prob( Chris@102: function, Chris@102: x_min, x_max, p, Chris@102: &result, Policy())) Chris@102: { Chris@102: return result; Chris@102: } Chris@102: // Special cases: Chris@102: if (p == 0) Chris@102: { Chris@102: return 0; Chris@102: } Chris@102: if (p == 1) Chris@102: { Chris@102: return 1; Chris@102: } Chris@102: Chris@102: RealType sin2hpip = sin(half_pi() * p); Chris@102: RealType sin2hpip2 = sin2hpip * sin2hpip; Chris@102: result = -x_min * sin2hpip2 + x_min + x_max * sin2hpip2; Chris@102: Chris@102: return result; Chris@102: } // quantile Chris@102: Chris@102: template Chris@102: inline RealType quantile(const complemented2_type, RealType>& c) Chris@102: { Chris@102: // Complement Quantile or Percent Point arcsine function. Chris@102: // Return the number of expected x for a given Chris@102: // complement of the probability q. Chris@102: BOOST_MATH_STD_USING // For ADL of std functions. Chris@102: Chris@102: using boost::math::constants::half_pi; Chris@102: static const char* function = "boost::math::quantile(arcsine_distribution<%1%> const&, %1%)"; Chris@102: Chris@102: // Error checks: Chris@102: RealType q = c.param; Chris@102: const arcsine_distribution& dist = c.dist; Chris@102: RealType result = 0; Chris@102: RealType x_min = dist.x_min(); Chris@102: RealType x_max = dist.x_max(); Chris@102: if (false == arcsine_detail::check_dist_and_prob( Chris@102: function, Chris@102: x_min, Chris@102: x_max, Chris@102: q, Chris@102: &result, Policy())) Chris@102: { Chris@102: return result; Chris@102: } Chris@102: // Special cases: Chris@102: if (q == 1) Chris@102: { Chris@102: return 0; Chris@102: } Chris@102: if (q == 0) Chris@102: { Chris@102: return 1; Chris@102: } Chris@102: // Naive RealType p = 1 - q; result = sin(half_pi() * p); loses accuracy, so use a cos alternative instead. Chris@102: //result = cos(half_pi() * q); // for arcsine(0,1) Chris@102: //result = result * result; Chris@102: // For generalized arcsine: Chris@102: RealType cos2hpip = cos(half_pi() * q); Chris@102: RealType cos2hpip2 = cos2hpip * cos2hpip; Chris@102: result = -x_min * cos2hpip2 + x_min + x_max * cos2hpip2; Chris@102: Chris@102: return result; Chris@102: } // Quantile Complement Chris@102: Chris@102: } // namespace math Chris@102: } // namespace boost Chris@102: Chris@102: // This include must be at the end, *after* the accessors Chris@102: // for this distribution have been defined, in order to Chris@102: // keep compilers that support two-phase lookup happy. Chris@102: #include Chris@102: Chris@102: #if defined (BOOST_MSVC) Chris@102: # pragma warning(pop) Chris@102: #endif Chris@102: Chris@102: #endif // BOOST_MATH_DIST_ARCSINE_HPP