diff --git a/doc/distributions/dist_reference.qbk b/doc/distributions/dist_reference.qbk index d8e3028ea6..de34e51d20 100644 --- a/doc/distributions/dist_reference.qbk +++ b/doc/distributions/dist_reference.qbk @@ -13,6 +13,7 @@ [include empirical_cdf.qbk] [include exponential.qbk] [include extreme_value.qbk] +[include generalized_extreme_value.qbk] [include fisher.qbk] [include gamma.qbk] [include geometric.qbk] diff --git a/doc/distributions/generalized_extreme_value.qbk b/doc/distributions/generalized_extreme_value.qbk new file mode 100644 index 0000000000..4b6980ae97 --- /dev/null +++ b/doc/distributions/generalized_extreme_value.qbk @@ -0,0 +1,116 @@ +[section:gev_dist Generalized Extreme Value Distribution] + +``#include `` + + template + class generalized_extreme_value_distribution; + + typedef generalized_extreme_value_distribution<> generalized_extreme_value; + + template + class generalized_extreme_value_distribution + { + public: + typedef RealType value_type; + + generalized_extreme_value_distribution(RealType location = 0, RealType scale = 1, RealType shape = 0); + + RealType location()const; + RealType scale()const; + RealType shape()const; + }; + +The [@https://en.wikipedia.org/wiki/Generalized_extreme_value_distribution generalized extreme value distribution] +is the family of limit distributions of the normalized maxima of independent, identically distributed +random variables. With location [mu], scale [sigma] and shape [xi], its CDF is + +[expression F(x) = e[super -t(x)], t(x) = (1 + [xi](x - [mu])/[sigma])[super -1/[xi]]] + +on the support where 1 + [xi](x - [mu])/[sigma] > 0, and its PDF is + +[expression f(x) = t(x)[super [xi]+1] e[super -t(x)] / [sigma]] + +The shape selects the type: + +* [xi] > 0, the Fr[eacute]chet type: heavy tailed to the right, and bounded below at [mu] - [sigma]/[xi]. +* [xi] = 0, the Gumbel type, where t(x) = e[super -(x-[mu])/[sigma]]: this is the __extreme_distrib. +* [xi] < 0, the reversed Weibull type: bounded above at [mu] - [sigma]/[xi]. + +Note that the sign convention for the shape varies: SciPy's `genextreme`, for example, takes c = -[xi]. + +[h4 Member Functions] + + generalized_extreme_value_distribution(RealType location = 0, RealType scale = 1, RealType shape = 0); + +Constructs a generalized extreme value distribution with the specified location, scale and shape. + +Requires that location and shape are finite and that scale is finite and > 0, otherwise calls __domain_error. + + RealType location()const; + RealType scale()const; + RealType shape()const; + +Return the parameters of the distribution. + +[h4 Non-member Accessors] + +All the [link math_toolkit.dist_ref.nmp usual non-member accessor functions] +that are generic to all distributions are supported: __usual_accessors. + +The domain of the random variable is \[-[infin], +[infin]\]: outside the support the PDF is zero +and the CDF is 0 or 1. `support` returns the bounded end for [xi] [ne] 0. + +The mean is infinite for [xi] [ge] 1 and the variance for [xi] [ge] 1/2, and both return +[infin] there. +The skewness is undefined for [xi] [ge] 1/3 and the kurtosis for [xi] [ge] 1/4, where they call __domain_error. The quantile at the unbounded ends calls __overflow_error, +as does the PDF at the upper end point when [xi] < -1, where the density is infinite. + +`logpdf`, `logcdf` and the log of the complement of the CDF are specialized for accuracy deep in the tails. + +[h4 Accuracy] + +Against references computed with mpmath at the double values of the arguments, the PDF, CDF, their +complements and logs, and the quantiles are within a few ulps, except where the function is itself +ill-conditioned: deep in the tails of the CDF, where e[super -t] magnifies the rounding of (x - [mu])/[sigma] +(up to about 30 ulps), and where the quantile is close to zero and [mu] + [sigma] q cancels. +The mean, variance, mode, median and entropy are within a few ulps. Skewness and kurtosis are accurate to +an absolute error of a few tens of epsilon for |[xi]| < 1/4, and of a few hundred for larger |[xi]|, +where they are third and fourth differences of log[Gamma](1 - k[xi]); both change sign +(the skewness near [xi] = -0.28, the excess kurtosis near [xi] = -0.15), so their relative error is unbounded there. + +[h4 Implementation] + +In the following table [mu] is the location, [sigma] the scale, [xi] the shape, z = (x - [mu])/[sigma], +/p/ the probability and /q = 1 - p/; g[sub k] = [Gamma](1 - k[xi]). + +[table +[[Function][Implementation Notes]] +[[pdf][exp(-(1 + [xi])L - t) / [sigma], where L = -log t = log1p([xi]z)/[xi], and t = e[super -L].]] +[[logpdf][-log [sigma] - (1 + [xi])L - t]] +[[cdf][exp(-t)]] +[[logcdf][-t]] +[[cdf complement][-expm1(-t)]] +[[log of the cdf complement][log(-expm1(-t)) for t < log 2, otherwise log1p(-e[super -t])]] +[[quantile][[mu] + [sigma] expm1(-[xi] log y)/[xi] with y = -log p, or [mu] + [sigma](y[super -[xi]] - 1)/[xi] when |[xi] log y| > 1.]] +[[quantile from the complement][As the quantile, with y = -log1p(-q).]] +[[mean][[mu] + [sigma] tgamma1pm1(-[xi])/[xi]]] +[[variance][([sigma]g[sub 1]/[xi])[super 2] expm1(A[sub 2]), with A[sub k] = log g[sub k] - k log g[sub 1]]] +[[skewness][sign([xi]) (g[sub 3] - 3g[sub 1]g[sub 2] + 2g[sub 1][super 3]) / (g[sub 2] - g[sub 1][super 2])[super 3/2], written in terms of the A[sub k] + so that the cancellation at small [xi] is exact: the combinations of the A[sub k] are summed from the Taylor series + log[Gamma](1 - z) = [gamma]z + [sum] [zeta](j) z[super j]/j, in which their low orders vanish identically, while the series converges quickly.]] +[[kurtosis excess][(g[sub 4] - 4g[sub 1]g[sub 3] + 6g[sub 1][super 2]g[sub 2] - 3g[sub 1][super 4]) / (g[sub 2] - g[sub 1][super 2])[super 2] - 3, as for the skewness.]] +[[mode][[mu] + [sigma]((1 + [xi])[super -[xi]] - 1)/[xi], or the upper end point [mu] - [sigma]/[xi] for [xi] [le] -1.]] +[[median][[mu] + [sigma]((log 2)[super -[xi]] - 1)/[xi]]] +[[entropy][log [sigma] + [gamma](1 + [xi]) + 1]] +] + +At [xi] = 0, and for |[xi]| below machine epsilon in the moments, the Gumbel formulas are used. + +[endsect] [/section:gev_dist Generalized Extreme Value] + +[/ generalized_extreme_value.qbk + Copyright 2026 Nick Thompson. + Distributed under the Boost Software License, Version 1.0. + (See accompanying file LICENSE_1_0.txt or copy at + http://www.boost.org/LICENSE_1_0.txt). +] diff --git a/include/boost/math/distributions.hpp b/include/boost/math/distributions.hpp index ca09602466..4ab64adb67 100644 --- a/include/boost/math/distributions.hpp +++ b/include/boost/math/distributions.hpp @@ -23,6 +23,7 @@ #include #include #include +#include #include #include #include diff --git a/include/boost/math/distributions/fwd.hpp b/include/boost/math/distributions/fwd.hpp index f166e57917..0be403a763 100644 --- a/include/boost/math/distributions/fwd.hpp +++ b/include/boost/math/distributions/fwd.hpp @@ -42,6 +42,9 @@ class extreme_value_distribution; BOOST_MATH_EXPORT template class fisher_f_distribution; +BOOST_MATH_EXPORT template +class generalized_extreme_value_distribution; + BOOST_MATH_EXPORT template class gamma_distribution; @@ -145,6 +148,7 @@ BOOST_MATH_NAMESPACE_END typedef BOOST_MATH_NAMESPACE::extreme_value_distribution extreme_value;\ typedef BOOST_MATH_NAMESPACE::fisher_f_distribution fisher_f;\ typedef BOOST_MATH_NAMESPACE::gamma_distribution gamma;\ + typedef BOOST_MATH_NAMESPACE::generalized_extreme_value_distribution generalized_extreme_value;\ typedef BOOST_MATH_NAMESPACE::geometric_distribution geometric;\ typedef BOOST_MATH_NAMESPACE::hypergeometric_distribution hypergeometric;\ typedef BOOST_MATH_NAMESPACE::kolmogorov_smirnov_distribution kolmogorov_smirnov;\ diff --git a/include/boost/math/distributions/generalized_extreme_value.hpp b/include/boost/math/distributions/generalized_extreme_value.hpp new file mode 100644 index 0000000000..722ac01718 --- /dev/null +++ b/include/boost/math/distributions/generalized_extreme_value.hpp @@ -0,0 +1,647 @@ +// Copyright Nick Thompson 2026. +// Use, modification and distribution are subject to the +// Boost Software License, Version 1.0. (See accompanying file +// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) + +#ifndef BOOST_MATH_DISTRIBUTIONS_GENERALIZED_EXTREME_VALUE_HPP +#define BOOST_MATH_DISTRIBUTIONS_GENERALIZED_EXTREME_VALUE_HPP + +// The generalized extreme value distribution, with location mu, scale sigma and shape xi as in +// https://en.wikipedia.org/wiki/Generalized_extreme_value_distribution: +// +// cdf(x) = exp(-t(x)), t(x) = (1 + xi (x - mu)/sigma)^(-1/xi), or exp(-(x - mu)/sigma) for xi = 0. +// +// xi > 0 is the Frechet type, bounded below at mu - sigma/xi; xi < 0 the reversed Weibull type, bounded +// above there; xi = 0 the Gumbel distribution, extreme_value_distribution. (SciPy's genextreme takes c = -xi.) +// log t = -log1p(xi z)/xi and the quantile's expm1(-xi log y)/xi tend to their Gumbel limits as xi -> 0 +// without cancellation. + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#ifndef BOOST_MATH_BUILD_MODULE +#include +#include +#include +#endif + +BOOST_MATH_NAMESPACE_BEGIN + +BOOST_MATH_EXPORT template > +class generalized_extreme_value_distribution +{ +public: + using value_type = RealType; + using policy_type = Policy; + + explicit generalized_extreme_value_distribution(RealType location = 0, RealType scale = 1, RealType shape = 0) + : m_location(location), m_scale(scale), m_shape(shape) + { + constexpr auto function = "boost::math::generalized_extreme_value_distribution<%1%>::generalized_extreme_value_distribution"; + RealType result; + detail::check_location(function, location, &result, Policy()); + detail::check_scale(function, scale, &result, Policy()); + detail::check_finite(function, shape, &result, Policy()); + } + + RealType location() const { return m_location; } + RealType scale() const { return m_scale; } + RealType shape() const { return m_shape; } + +private: + RealType m_location; + RealType m_scale; + RealType m_shape; +}; + +BOOST_MATH_EXPORT using generalized_extreme_value = generalized_extreme_value_distribution; + +#ifdef __cpp_deduction_guides +BOOST_MATH_EXPORT template +generalized_extreme_value_distribution(RealType)->generalized_extreme_value_distribution::type>; +BOOST_MATH_EXPORT template +generalized_extreme_value_distribution(RealType, RealType)->generalized_extreme_value_distribution::type>; +BOOST_MATH_EXPORT template +generalized_extreme_value_distribution(RealType, RealType, RealType)->generalized_extreme_value_distribution::type>; +#endif + +namespace detail { + +template +inline bool check_gev(const char* function, const generalized_extreme_value_distribution& dist, RealType* result) +{ + return check_location(function, dist.location(), result, Policy()) + && check_scale(function, dist.scale(), result, Policy()) + && check_finite(function, dist.shape(), result, Policy()); +} + +// Where x lies relative to the support, and inside it L = -log t(x) = log1p(xi z)/xi, so that cdf(x) = exp(-exp(-L)). +enum class gev_region { below, inside, endpoint, above }; + +// The support is decided against the end point mu - sigma/xi that support() reports, before x is standardized, +// whose rounding could otherwise carry a point outside onto it. Inside, 1 + xi z can still round to 0 or below +// next to the end point, where the end point's limit is the answer to working precision. +template +gev_region gev_log_t(RealType x, RealType mu, RealType sigma, RealType xi, RealType* L) +{ + BOOST_MATH_STD_USING + RealType z = (x - mu) / sigma; + if (xi == 0) + { + *L = z; + return gev_region::inside; + } + RealType endpoint = mu - sigma / xi; + if (x == endpoint) + { + return gev_region::endpoint; + } + if ((xi > 0) == (x < endpoint)) + { + // Beyond the lower end point for xi > 0, the upper for xi < 0: + return xi > 0 ? gev_region::below : gev_region::above; + } + RealType s = xi * z; + if (s <= -1) + { + return gev_region::endpoint; + } + // For tiny xi z, xi z may have lost bits to underflow; log1p(s)/xi = z (1 - s/2 + ...): + *L = abs(s) < tools::epsilon() ? z * (1 - s / 2) : RealType(BOOST_MATH_NAMESPACE::log1p(s) / xi); + return gev_region::inside; +} + +// (y^(-xi) - 1)/xi, the standardized quantile at t = y; -log y for xi = 0. expm1(-xi log y) carries +// log's absolute error times xi, so where |xi log y| > 1 pow, which then has nothing to cancel, is better. +template +RealType gev_standard_quantile(RealType y, RealType log_y, RealType xi, const Policy& pol) +{ + BOOST_MATH_STD_USING + RealType a = -xi * log_y; + if (xi == 0 || abs(a) < tools::epsilon()) + { + return -log_y * (1 + a / 2); + } + if (abs(a) > 1) + { + return (pow(y, -xi) - 1) / xi; + } + return BOOST_MATH_NAMESPACE::expm1(a, pol) / xi; +} + +// expm1(x) - x, without the cancellation for small x. +template +RealType gev_expm1mx(RealType x, const Policy& pol) +{ + BOOST_MATH_STD_USING + if (abs(x) > RealType(0.5)) + { + return BOOST_MATH_NAMESPACE::expm1(x, pol) - x; + } + RealType term = x * x / 2; + RealType sum = term; + for (unsigned n = 3; abs(term) > tools::epsilon() * abs(sum); ++n) + { + term *= x / n; + sum += term; + } + return sum; +} + +// sum_k c[k] A_k for k = 2, 3, 4, where A_k = lgamma(1 - k xi) - k lgamma(1 - xi) is the log of +// E[W^k]/E[W]^k for W = t^(-xi), the moments of which are Gamma(1 - k xi). The combinations the +// moments need cancel at leading order(s) in xi, so for small xi this sums the Taylor series +// lgamma(1 - z) = euler z + sum_{j >= 2} zeta(j) z^j / j term by term, in which they cancel exactly: +// sum_k c[k] A_k = sum_{j >= 2} zeta(j)/j sum_k c[k] (k^j - k) xi^j. +// The series converges for |k xi| < 1 with ratio |k xi|; it is used to |k xi| = 3/4 for the largest k, +// beyond which the lgamma differences lose at most a digit or so. +template +RealType gev_log_moment_combination(const int (&c)[3], RealType xi, const Policy& pol) +{ + BOOST_MATH_STD_USING + int k_max = c[2] != 0 ? 4 : (c[1] != 0 ? 3 : 2); + if (k_max * abs(xi) > RealType(0.75)) + { + RealType l1 = BOOST_MATH_NAMESPACE::lgamma(1 - xi, pol); + RealType sum = 0; + for (int k = 2; k <= 4; ++k) + { + if (c[k - 2] != 0) + { + sum += c[k - 2] * (BOOST_MATH_NAMESPACE::lgamma(1 - k * xi, pol) - k * l1); + } + } + return sum; + } + RealType power[3] = { 2, 3, 4 }; // k^j + RealType xi_j = xi; // xi^j + RealType sum = 0; + const std::uintmax_t max_terms = policies::get_max_series_iterations(); + for (std::uintmax_t j = 2; ; ++j) + { + RealType weight = 0; + for (int k = 2; k <= 4; ++k) + { + power[k - 2] *= k; + weight += c[k - 2] * (power[k - 2] - k); + } + xi_j *= xi; + RealType term = BOOST_MATH_NAMESPACE::zeta(RealType(j), pol) / j * weight * xi_j; + sum += term; + // The weights of the leading orders are zero, so only stop on a nonzero term: + if (weight != 0 && abs(term) <= tools::epsilon() * abs(sum)) + { + return sum; + } + if (j > max_terms) + { + return policies::raise_evaluation_error("boost::math::generalized_extreme_value moments", "Series did not converge for shape %1%.", xi, pol); + } + } +} + +// The moments of xi this small are their Gumbel limits to working precision. +template +bool gev_is_gumbel(RealType xi) +{ + BOOST_MATH_STD_USING + return abs(xi) < tools::epsilon(); +} + +// The mean for xi >= 1 and the variance for xi >= 1/2, which diverge to +infinity (max_value for a type without one). +template +RealType gev_infinity() +{ + return std::numeric_limits::has_infinity ? std::numeric_limits::infinity() : tools::max_value(); +} + +template +RealType gev_undefined_moment(const char* function, const char* message, RealType xi) +{ + return policies::raise_domain_error(function, message, xi, Policy()); +} + +} // namespace detail + +BOOST_MATH_EXPORT template +inline std::pair range(const generalized_extreme_value_distribution& /*dist*/) +{ // Range of permissible values for random variable x; the pdf is zero outside the support. + using BOOST_MATH_NAMESPACE::tools::max_value; + return std::pair( + std::numeric_limits::has_infinity ? -std::numeric_limits::infinity() : -max_value(), + std::numeric_limits::has_infinity ? std::numeric_limits::infinity() : max_value()); +} + +BOOST_MATH_EXPORT template +inline std::pair support(const generalized_extreme_value_distribution& dist) +{ // Range where the cdf rises from 0 to 1: bounded at mu - sigma/xi below for xi > 0, above for xi < 0. + using BOOST_MATH_NAMESPACE::tools::max_value; + RealType xi = dist.shape(); + RealType endpoint = xi == 0 ? RealType(0) : RealType(dist.location() - dist.scale() / xi); + if (xi > 0) + { + return std::pair(endpoint, max_value()); + } + if (xi < 0) + { + return std::pair(-max_value(), endpoint); + } + return std::pair(-max_value(), max_value()); +} + +BOOST_MATH_EXPORT template +inline RealType logpdf(const generalized_extreme_value_distribution& dist, const RealType& x) +{ + BOOST_MATH_STD_USING + constexpr auto function = "boost::math::logpdf(const generalized_extreme_value_distribution<%1%>&, %1%)"; + RealType result = 0; + if (!detail::check_gev(function, dist, &result) || !detail::check_x_not_NaN(function, x, &result, Policy())) + { + return result; + } + RealType xi = dist.shape(); + if ((BOOST_MATH_NAMESPACE::isinf)(x)) + { + return -std::numeric_limits::infinity(); + } + RealType L; + switch (detail::gev_log_t(x, dist.location(), dist.scale(), xi, &L)) + { + case detail::gev_region::inside: + break; + case detail::gev_region::endpoint: + // The bounded end, where t = 0 for xi < 0 and the density is t^(1 + xi)/sigma; t = infinity for xi > 0. + if (xi == -1) + { + return -log(dist.scale()); + } + if (xi < -1) + { + return policies::raise_overflow_error(function, nullptr, Policy()); + } + return -std::numeric_limits::infinity(); + default: + return -std::numeric_limits::infinity(); + } + // log pdf = -log sigma + (1 + xi) log t - t, with log t = -L: + RealType t = exp(-L); + if ((BOOST_MATH_NAMESPACE::isinf)(t)) + { + return -std::numeric_limits::infinity(); + } + return -log(dist.scale()) - (1 + xi) * L - t; +} + +BOOST_MATH_EXPORT template +inline RealType pdf(const generalized_extreme_value_distribution& dist, const RealType& x) +{ + BOOST_MATH_STD_USING + constexpr auto function = "boost::math::pdf(const generalized_extreme_value_distribution<%1%>&, %1%)"; + RealType result = 0; + if (!detail::check_gev(function, dist, &result) || !detail::check_x_not_NaN(function, x, &result, Policy())) + { + return result; + } + RealType xi = dist.shape(); + if ((BOOST_MATH_NAMESPACE::isinf)(x)) + { + return 0; + } + RealType L; + switch (detail::gev_log_t(x, dist.location(), dist.scale(), xi, &L)) + { + case detail::gev_region::inside: + break; + case detail::gev_region::endpoint: + if (xi == -1) + { + return 1 / dist.scale(); + } + if (xi < -1) + { + return policies::raise_overflow_error(function, nullptr, Policy()); + } + return 0; + default: + return 0; + } + RealType t = exp(-L); + if ((BOOST_MATH_NAMESPACE::isinf)(t)) + { + return 0; + } + // t^(1 + xi) exp(-t) / sigma: + return exp(-(1 + xi) * L - t) / dist.scale(); +} + +namespace detail { + +// cdf, its complement, and their logs; region and t decide them all. +template +RealType gev_cdf_imp(const generalized_extreme_value_distribution& dist, const RealType& x, bool complement, bool take_log, const char* function) +{ + BOOST_MATH_STD_USING + RealType result = 0; + if (!check_gev(function, dist, &result) || !check_x_not_NaN(function, x, &result, Policy())) + { + return result; + } + const RealType zero = take_log ? RealType(-std::numeric_limits::infinity()) : RealType(0); + const RealType one = take_log ? RealType(0) : RealType(1); + if ((BOOST_MATH_NAMESPACE::isinf)(x)) + { + return (x < 0) != complement ? zero : one; + } + RealType xi = dist.shape(); + RealType L; + switch (gev_log_t(x, dist.location(), dist.scale(), xi, &L)) + { + case gev_region::inside: + break; + case gev_region::endpoint: + return (xi > 0) != complement ? zero : one; + case gev_region::below: + return complement ? one : zero; + default: + return complement ? zero : one; + } + RealType t = exp(-L); + if (!complement) + { + return take_log ? RealType(-t) : RealType(exp(-t)); + } + RealType q = -BOOST_MATH_NAMESPACE::expm1(-t, Policy()); + if (!take_log) + { + return q; + } + // log(1 - exp(-t)): log t - t/2 = -L - t/2 once t is below epsilon (t may have underflowed), through q + // while it is small, else through log1p: + if (t < tools::epsilon()) + { + return -L - t / 2; + } + return t () ? RealType(log(q)) : RealType(BOOST_MATH_NAMESPACE::log1p(-exp(-t), Policy())); +} + +template +RealType gev_quantile_imp(const generalized_extreme_value_distribution& dist, const RealType& p, bool complement, const char* function) +{ + BOOST_MATH_STD_USING + RealType result = 0; + if (!check_gev(function, dist, &result) || !check_probability(function, p, &result, Policy())) + { + return result; + } + RealType xi = dist.shape(); + RealType mu = dist.location(); + RealType sigma = dist.scale(); + // P = cdf(x): + bool lower_end = complement ? p == 1 : p == 0; + bool upper_end = complement ? p == 0 : p == 1; + if (lower_end) + { + return xi > 0 ? RealType(mu - sigma / xi) : RealType(-policies::raise_overflow_error(function, nullptr, Policy())); + } + if (upper_end) + { + return xi < 0 ? RealType(mu - sigma / xi) : policies::raise_overflow_error(function, nullptr, Policy()); + } + // t = -log P: + RealType y = complement ? RealType(-BOOST_MATH_NAMESPACE::log1p(-p, Policy())) : RealType(-log(p)); + return mu + sigma * gev_standard_quantile(y, RealType(log(y)), xi, Policy()); +} + +} // namespace detail + +BOOST_MATH_EXPORT template +inline RealType cdf(const generalized_extreme_value_distribution& dist, const RealType& x) +{ + return detail::gev_cdf_imp(dist, x, false, false, "boost::math::cdf(const generalized_extreme_value_distribution<%1%>&, %1%)"); +} + +BOOST_MATH_EXPORT template +inline RealType cdf(const complemented2_type, RealType>& c) +{ + return detail::gev_cdf_imp(c.dist, c.param, true, false, "boost::math::cdf(const generalized_extreme_value_distribution<%1%>&, %1%)"); +} + +BOOST_MATH_EXPORT template +inline RealType logcdf(const generalized_extreme_value_distribution& dist, const RealType& x) +{ + return detail::gev_cdf_imp(dist, x, false, true, "boost::math::logcdf(const generalized_extreme_value_distribution<%1%>&, %1%)"); +} + +BOOST_MATH_EXPORT template +inline RealType logcdf(const complemented2_type, RealType>& c) +{ + return detail::gev_cdf_imp(c.dist, c.param, true, true, "boost::math::logcdf(const generalized_extreme_value_distribution<%1%>&, %1%)"); +} + +BOOST_MATH_EXPORT template +inline RealType quantile(const generalized_extreme_value_distribution& dist, const RealType& p) +{ + return detail::gev_quantile_imp(dist, p, false, "boost::math::quantile(const generalized_extreme_value_distribution<%1%>&, %1%)"); +} + +BOOST_MATH_EXPORT template +inline RealType quantile(const complemented2_type, RealType>& c) +{ + return detail::gev_quantile_imp(c.dist, c.param, true, "boost::math::quantile(const generalized_extreme_value_distribution<%1%>&, %1%)"); +} + +BOOST_MATH_EXPORT template +inline RealType mean(const generalized_extreme_value_distribution& dist) +{ + constexpr auto function = "boost::math::mean(const generalized_extreme_value_distribution<%1%>&)"; + RealType result = 0; + if (!detail::check_gev(function, dist, &result)) + { + return result; + } + RealType xi = dist.shape(); + if (xi >= 1) + { + // The right tail decays like x^(-1/xi), so the mean is +infinity: + return detail::gev_infinity(); + } + if (detail::gev_is_gumbel(xi)) + { + return mean(extreme_value_distribution(dist.location(), dist.scale())); + } + // mu + sigma (Gamma(1 - xi) - 1)/xi: + return dist.location() + dist.scale() * BOOST_MATH_NAMESPACE::tgamma1pm1(-xi, Policy()) / xi; +} + +BOOST_MATH_EXPORT template +inline RealType variance(const generalized_extreme_value_distribution& dist) +{ + constexpr auto function = "boost::math::variance(const generalized_extreme_value_distribution<%1%>&)"; + RealType result = 0; + if (!detail::check_gev(function, dist, &result)) + { + return result; + } + RealType xi = dist.shape(); + if (2 * xi >= 1) + { + // The second moment diverges, so the variance is +infinity: + return detail::gev_infinity(); + } + RealType sigma = dist.scale(); + if (detail::gev_is_gumbel(xi)) + { + return sigma * sigma * constants::pi_sqr() / 6; + } + // sigma^2 (Gamma(1 - 2 xi) - Gamma(1 - xi)^2)/xi^2 = (sigma Gamma(1 - xi)/xi)^2 expm1(A_2): + const int a2[3] = { 1, 0, 0 }; + RealType g1 = BOOST_MATH_NAMESPACE::tgamma(1 - xi, Policy()); + RealType s = sigma * g1 / xi; + return s * s * BOOST_MATH_NAMESPACE::expm1(detail::gev_log_moment_combination(a2, xi, Policy()), Policy()); +} + +BOOST_MATH_EXPORT template +inline RealType standard_deviation(const generalized_extreme_value_distribution& dist) +{ + BOOST_MATH_STD_USING + return sqrt(variance(dist)); +} + +BOOST_MATH_EXPORT template +inline RealType skewness(const generalized_extreme_value_distribution& dist) +{ + BOOST_MATH_STD_USING + constexpr auto function = "boost::math::skewness(const generalized_extreme_value_distribution<%1%>&)"; + RealType result = 0; + if (!detail::check_gev(function, dist, &result)) + { + return result; + } + RealType xi = dist.shape(); + if (3 * xi >= 1) + { + return detail::gev_undefined_moment(function, "The skewness is undefined for shape >= 1/3, but got %1%.", xi); + } + if (detail::gev_is_gumbel(xi)) + { + return skewness(extreme_value_distribution()); + } + // sign(xi) (g3 - 3 g1 g2 + 2 g1^3)/(g2 - g1^2)^(3/2) with g_k = Gamma(1 - k xi); over g1^3 the numerator + // is expm1(A_3) - 3 expm1(A_2) = (A_3 - 3 A_2) + expm1mx(A_3) - 3 expm1mx(A_2), which is O(xi^3): + const int a2[3] = { 1, 0, 0 }; + const int a3[3] = { 0, 1, 0 }; + const int b[3] = { -3, 1, 0 }; + RealType A2 = detail::gev_log_moment_combination(a2, xi, Policy()); + RealType A3 = detail::gev_log_moment_combination(a3, xi, Policy()); + RealType numerator = detail::gev_log_moment_combination(b, xi, Policy()) + detail::gev_expm1mx(A3, Policy()) - 3 * detail::gev_expm1mx(A2, Policy()); + RealType v = BOOST_MATH_NAMESPACE::expm1(A2, Policy()); + RealType skew = numerator / (v * sqrt(v)); + return xi > 0 ? skew : RealType(-skew); +} + +BOOST_MATH_EXPORT template +inline RealType kurtosis_excess(const generalized_extreme_value_distribution& dist) +{ + constexpr auto function = "boost::math::kurtosis_excess(const generalized_extreme_value_distribution<%1%>&)"; + RealType result = 0; + if (!detail::check_gev(function, dist, &result)) + { + return result; + } + RealType xi = dist.shape(); + if (4 * xi >= 1) + { + return detail::gev_undefined_moment(function, "The kurtosis is undefined for shape >= 1/4, but got %1%.", xi); + } + if (detail::gev_is_gumbel(xi)) + { + return kurtosis_excess(extreme_value_distribution()); + } + // (g4 - 4 g1 g3 + 6 g1^2 g2 - 3 g1^4)/(g2 - g1^2)^2 - 3; over g1^4 the numerator is + // expm1(A_4) - 4 expm1(A_3) + 6 expm1(A_2), whose part linear in the A_k is O(xi^4): + const int a2[3] = { 1, 0, 0 }; + const int a3[3] = { 0, 1, 0 }; + const int a4[3] = { 0, 0, 1 }; + const int c[3] = { 6, -4, 1 }; + RealType A2 = detail::gev_log_moment_combination(a2, xi, Policy()); + RealType A3 = detail::gev_log_moment_combination(a3, xi, Policy()); + RealType A4 = detail::gev_log_moment_combination(a4, xi, Policy()); + RealType numerator = detail::gev_log_moment_combination(c, xi, Policy()) + + detail::gev_expm1mx(A4, Policy()) - 4 * detail::gev_expm1mx(A3, Policy()) + 6 * detail::gev_expm1mx(A2, Policy()); + RealType v = BOOST_MATH_NAMESPACE::expm1(A2, Policy()); + return numerator / (v * v) - 3; +} + +BOOST_MATH_EXPORT template +inline RealType kurtosis(const generalized_extreme_value_distribution& dist) +{ + return kurtosis_excess(dist) + 3; +} + +BOOST_MATH_EXPORT template +inline RealType mode(const generalized_extreme_value_distribution& dist) +{ + BOOST_MATH_STD_USING + constexpr auto function = "boost::math::mode(const generalized_extreme_value_distribution<%1%>&)"; + RealType result = 0; + if (!detail::check_gev(function, dist, &result)) + { + return result; + } + RealType xi = dist.shape(); + if (xi <= -1) + { + // The density increases to the upper endpoint: + return dist.location() - dist.scale() / xi; + } + // mu + sigma ((1 + xi)^(-xi) - 1)/xi, where t = 1 + xi: + return dist.location() + dist.scale() * detail::gev_standard_quantile(RealType(1 + xi), RealType(BOOST_MATH_NAMESPACE::log1p(xi, Policy())), xi, Policy()); +} + +BOOST_MATH_EXPORT template +inline RealType median(const generalized_extreme_value_distribution& dist) +{ + constexpr auto function = "boost::math::median(const generalized_extreme_value_distribution<%1%>&)"; + RealType result = 0; + if (!detail::check_gev(function, dist, &result)) + { + return result; + } + // t = log 2: + return dist.location() + dist.scale() * detail::gev_standard_quantile(constants::ln_two(), constants::ln_ln_two(), dist.shape(), Policy()); +} + +BOOST_MATH_EXPORT template +inline RealType entropy(const generalized_extreme_value_distribution& dist) +{ + BOOST_MATH_STD_USING + constexpr auto function = "boost::math::entropy(const generalized_extreme_value_distribution<%1%>&)"; + RealType result = 0; + if (!detail::check_gev(function, dist, &result)) + { + return result; + } + return log(dist.scale()) + constants::euler() * (1 + dist.shape()) + 1; +} + +BOOST_MATH_NAMESPACE_END + +// This include must be at the end, *after* the accessors +// for this distribution have been defined, in order to +// keep compilers that support two-phase lookup happy. +#include + +#endif // BOOST_MATH_DISTRIBUTIONS_GENERALIZED_EXTREME_VALUE_HPP diff --git a/nanobind/src/distributions.hpp b/nanobind/src/distributions.hpp index 6f6e5f65fe..44852970c1 100644 --- a/nanobind/src/distributions.hpp +++ b/nanobind/src/distributions.hpp @@ -14,7 +14,10 @@ #include +#include +#include #include +#include #include "vectorize.hpp" @@ -70,6 +73,35 @@ nb::class_& add_dist_methods(nb::class_& cls) { return cls; } +// A parameter: the constructor's argument name, and its value read back from the distribution. +template +using dist_param = std::pair; + +#define BMPY_PARAM(name, accessor) \ + bmpy::dist_param{name, +[](const D& d) { return nb::cast(d.accessor()); }} + +// Exposes the parameters as read-only properties, under the constructor's argument names, and a repr that +// constructs the same distribution: normal_distribution(mean=0.0, sd=1.0). A parameter named like a method +// (the mean of a normal, the mode of a triangular) stays the method, which returns the same value. +template +nb::class_& add_dist_params(nb::class_& cls, std::initializer_list> params) { + std::vector> ps(params); + for (auto const& [name, get] : ps) { + if (!nb::hasattr(cls, name)) { + cls.def_prop_ro(name, get); + } + } + std::string type = nb::cast(cls.attr("__name__")); + cls.def("__repr__", [type, ps](const D& d) { + std::string s = type + "("; + for (std::size_t i = 0; i < ps.size(); ++i) { + s += (i ? ", " : "") + std::string(ps[i].first) + "=" + nb::cast(nb::repr(ps[i].second(d))); + } + return s + ")"; + }); + return cls; +} + } // namespace bmpy #endif diff --git a/nanobind/src/distributions_continuous.cpp b/nanobind/src/distributions_continuous.cpp index 87757c0927..4c9a6c51fb 100644 --- a/nanobind/src/distributions_continuous.cpp +++ b/nanobind/src/distributions_continuous.cpp @@ -10,6 +10,7 @@ #include #include #include +#include #include #include #include @@ -40,48 +41,63 @@ void bmpy::bind_distributions_continuous(nb::module_& m) { nb::class_ c(m, "arcsine_distribution"); c.def(nb::init(), arg("x_min") = 0.0, arg("x_max") = 1.0); add_dist_methods(c); + bmpy::add_dist_params(c, {BMPY_PARAM("x_min", x_min), BMPY_PARAM("x_max", x_max)}); } { using D = bm::beta_distribution; nb::class_ c(m, "beta_distribution"); c.def(nb::init(), arg("alpha") = 1.0, arg("beta") = 1.0); add_dist_methods(c); + bmpy::add_dist_params(c, {BMPY_PARAM("alpha", alpha), BMPY_PARAM("beta", beta)}); } { using D = bm::cauchy_distribution; nb::class_ c(m, "cauchy_distribution"); c.def(nb::init(), arg("location") = 0.0, arg("scale") = 1.0); add_dist_methods(c); + bmpy::add_dist_params(c, {BMPY_PARAM("location", location), BMPY_PARAM("scale", scale)}); } { using D = bm::chi_squared_distribution; nb::class_ c(m, "chi_squared_distribution"); c.def(nb::init(), arg("df") = 1.0); add_dist_methods(c); + bmpy::add_dist_params(c, {BMPY_PARAM("df", degrees_of_freedom)}); } { using D = bm::exponential_distribution; nb::class_ c(m, "exponential_distribution"); c.def(nb::init(), arg("lambda_") = 1.0); add_dist_methods(c); + bmpy::add_dist_params(c, {BMPY_PARAM("lambda_", lambda)}); } { using D = bm::extreme_value_distribution; nb::class_ c(m, "extreme_value_distribution"); c.def(nb::init(), arg("location") = 0.0, arg("scale") = 1.0); add_dist_methods(c); + bmpy::add_dist_params(c, {BMPY_PARAM("location", location), BMPY_PARAM("scale", scale)}); + } + { + using D = bm::generalized_extreme_value_distribution; + nb::class_ c(m, "generalized_extreme_value_distribution"); + c.def(nb::init(), arg("location") = 0.0, arg("scale") = 1.0, arg("shape") = 0.0); + add_dist_methods(c); + bmpy::add_dist_params(c, {BMPY_PARAM("location", location), BMPY_PARAM("scale", scale), BMPY_PARAM("shape", shape)}); } { using D = bm::fisher_f_distribution; nb::class_ c(m, "fisher_f_distribution"); c.def(nb::init(), arg("df1"), arg("df2")); add_dist_methods(c); + bmpy::add_dist_params(c, {BMPY_PARAM("df1", degrees_of_freedom1), BMPY_PARAM("df2", degrees_of_freedom2)}); } { using D = bm::gamma_distribution; nb::class_ c(m, "gamma_distribution"); c.def(nb::init(), arg("shape"), arg("scale") = 1.0); add_dist_methods(c); + bmpy::add_dist_params(c, {BMPY_PARAM("shape", shape), BMPY_PARAM("scale", scale)}); } { using D = bm::inverse_chi_squared_distribution; @@ -89,83 +105,97 @@ void bmpy::bind_distributions_continuous(nb::module_& m) { c.def(nb::init(), arg("df") = 1.0); c.def(nb::init(), arg("df"), arg("scale")); add_dist_methods(c); + bmpy::add_dist_params(c, {BMPY_PARAM("df", degrees_of_freedom), BMPY_PARAM("scale", scale)}); } { using D = bm::inverse_gamma_distribution; nb::class_ c(m, "inverse_gamma_distribution"); c.def(nb::init(), arg("shape") = 1.0, arg("scale") = 1.0); add_dist_methods(c); + bmpy::add_dist_params(c, {BMPY_PARAM("shape", shape), BMPY_PARAM("scale", scale)}); } { using D = bm::inverse_gaussian_distribution; nb::class_ c(m, "inverse_gaussian_distribution"); c.def(nb::init(), arg("mean") = 1.0, arg("scale") = 1.0); add_dist_methods(c); + bmpy::add_dist_params(c, {BMPY_PARAM("mean", mean), BMPY_PARAM("scale", scale)}); } { using D = bm::laplace_distribution; nb::class_ c(m, "laplace_distribution"); c.def(nb::init(), arg("location") = 0.0, arg("scale") = 1.0); add_dist_methods(c); + bmpy::add_dist_params(c, {BMPY_PARAM("location", location), BMPY_PARAM("scale", scale)}); } { using D = bm::logistic_distribution; nb::class_ c(m, "logistic_distribution"); c.def(nb::init(), arg("location") = 0.0, arg("scale") = 1.0); add_dist_methods(c); + bmpy::add_dist_params(c, {BMPY_PARAM("location", location), BMPY_PARAM("scale", scale)}); } { using D = bm::lognormal_distribution; nb::class_ c(m, "lognormal_distribution"); c.def(nb::init(), arg("location") = 0.0, arg("scale") = 1.0); add_dist_methods(c); + bmpy::add_dist_params(c, {BMPY_PARAM("location", location), BMPY_PARAM("scale", scale)}); } { using D = bm::normal_distribution; nb::class_ c(m, "normal_distribution"); c.def(nb::init(), arg("mean") = 0.0, arg("sd") = 1.0); add_dist_methods(c); + bmpy::add_dist_params(c, {BMPY_PARAM("mean", mean), BMPY_PARAM("sd", standard_deviation)}); } { using D = bm::pareto_distribution; nb::class_ c(m, "pareto_distribution"); c.def(nb::init(), arg("scale") = 1.0, arg("shape") = 1.0); add_dist_methods(c); + bmpy::add_dist_params(c, {BMPY_PARAM("scale", scale), BMPY_PARAM("shape", shape)}); } { using D = bm::rayleigh_distribution; nb::class_ c(m, "rayleigh_distribution"); c.def(nb::init(), arg("sigma") = 1.0); add_dist_methods(c); + bmpy::add_dist_params(c, {BMPY_PARAM("sigma", sigma)}); } { using D = bm::skew_normal_distribution; nb::class_ c(m, "skew_normal_distribution"); c.def(nb::init(), arg("location") = 0.0, arg("scale") = 1.0, arg("shape") = 0.0); add_dist_methods(c); + bmpy::add_dist_params(c, {BMPY_PARAM("location", location), BMPY_PARAM("scale", scale), BMPY_PARAM("shape", shape)}); } { using D = bm::students_t_distribution; nb::class_ c(m, "students_t_distribution"); c.def(nb::init(), arg("df") = 1.0); add_dist_methods(c); + bmpy::add_dist_params(c, {BMPY_PARAM("df", degrees_of_freedom)}); } { using D = bm::triangular_distribution; nb::class_ c(m, "triangular_distribution"); c.def(nb::init(), arg("lower") = -1.0, arg("mode") = 0.0, arg("upper") = 1.0); add_dist_methods(c); + bmpy::add_dist_params(c, {BMPY_PARAM("lower", lower), BMPY_PARAM("mode", mode), BMPY_PARAM("upper", upper)}); } { using D = bm::uniform_distribution; nb::class_ c(m, "uniform_distribution"); c.def(nb::init(), arg("lower") = 0.0, arg("upper") = 1.0); add_dist_methods(c); + bmpy::add_dist_params(c, {BMPY_PARAM("lower", lower), BMPY_PARAM("upper", upper)}); } { using D = bm::weibull_distribution; nb::class_ c(m, "weibull_distribution"); c.def(nb::init(), arg("shape"), arg("scale") = 1.0); add_dist_methods(c); + bmpy::add_dist_params(c, {BMPY_PARAM("shape", shape), BMPY_PARAM("scale", scale)}); } } diff --git a/nanobind/src/distributions_discrete.cpp b/nanobind/src/distributions_discrete.cpp index 6aa7ad61ec..bb8cf5354a 100644 --- a/nanobind/src/distributions_discrete.cpp +++ b/nanobind/src/distributions_discrete.cpp @@ -28,35 +28,41 @@ void bmpy::bind_distributions_discrete(nb::module_& m) { nb::class_ c(m, "bernoulli_distribution"); c.def(nb::init(), arg("p") = 0.5); add_dist_methods(c); + bmpy::add_dist_params(c, {BMPY_PARAM("p", success_fraction)}); } { using D = bm::binomial_distribution; nb::class_ c(m, "binomial_distribution"); c.def(nb::init(), arg("n") = 1.0, arg("p") = 0.5); add_dist_methods(c); + bmpy::add_dist_params(c, {BMPY_PARAM("n", trials), BMPY_PARAM("p", success_fraction)}); } { using D = bm::geometric_distribution; nb::class_ c(m, "geometric_distribution"); c.def(nb::init(), arg("p")); add_dist_methods(c); + bmpy::add_dist_params(c, {BMPY_PARAM("p", success_fraction)}); } { using D = bm::hypergeometric_distribution; nb::class_ c(m, "hypergeometric_distribution"); c.def(nb::init(), arg("r"), arg("n"), arg("N")); add_dist_methods(c); + bmpy::add_dist_params(c, {BMPY_PARAM("r", defective), BMPY_PARAM("n", sample_count), BMPY_PARAM("N", total)}); } { using D = bm::negative_binomial_distribution; nb::class_ c(m, "negative_binomial_distribution"); c.def(nb::init(), arg("r"), arg("p")); add_dist_methods(c); + bmpy::add_dist_params(c, {BMPY_PARAM("r", successes), BMPY_PARAM("p", success_fraction)}); } { using D = bm::poisson_distribution; nb::class_ c(m, "poisson_distribution"); c.def(nb::init(), arg("mean") = 1.0); add_dist_methods(c); + bmpy::add_dist_params(c, {BMPY_PARAM("mean", mean)}); } } diff --git a/nanobind/src/distributions_heavy.cpp b/nanobind/src/distributions_heavy.cpp index 45764e7dd1..3687b8c569 100644 --- a/nanobind/src/distributions_heavy.cpp +++ b/nanobind/src/distributions_heavy.cpp @@ -31,36 +31,42 @@ void bmpy::bind_distributions_heavy(nb::module_& m) { nb::class_ c(m, "holtsmark_distribution"); c.def(nb::init(), arg("location") = 0.0, arg("scale") = 1.0); add_dist_methods(c); + bmpy::add_dist_params(c, {BMPY_PARAM("location", location), BMPY_PARAM("scale", scale)}); } { using D = bm::landau_distribution; nb::class_ c(m, "landau_distribution"); c.def(nb::init(), arg("location") = 0.0, arg("scale") = 1.0); add_dist_methods(c); + bmpy::add_dist_params(c, {BMPY_PARAM("location", location), BMPY_PARAM("scale", scale)}); } { using D = bm::mapairy_distribution; nb::class_ c(m, "mapairy_distribution"); c.def(nb::init(), arg("location") = 0.0, arg("scale") = 1.0); add_dist_methods(c); + bmpy::add_dist_params(c, {BMPY_PARAM("location", location), BMPY_PARAM("scale", scale)}); } { using D = bm::saspoint5_distribution; nb::class_ c(m, "saspoint5_distribution"); c.def(nb::init(), arg("location") = 0.0, arg("scale") = 1.0); add_dist_methods(c); + bmpy::add_dist_params(c, {BMPY_PARAM("location", location), BMPY_PARAM("scale", scale)}); } { using D = bm::von_mises_distribution; nb::class_ c(m, "von_mises_distribution"); c.def(nb::init(), arg("mean") = 0.0, arg("concentration") = 1.0); add_dist_methods(c); + bmpy::add_dist_params(c, {BMPY_PARAM("mean", mean), BMPY_PARAM("concentration", concentration)}); } { using D = bm::kolmogorov_smirnov_distribution; nb::class_ c(m, "kolmogorov_smirnov_distribution"); c.def(nb::init(), arg("n")); add_dist_methods(c); + bmpy::add_dist_params(c, {BMPY_PARAM("n", number_of_observations)}); } { using D = bm::hyperexponential_distribution; @@ -71,5 +77,6 @@ void bmpy::bind_distributions_heavy(nb::module_& m) { [](D* self, const std::vector& probs, const std::vector& rates) { new (self) D(probs, rates); }, arg("probabilities"), arg("rates")); add_dist_methods(c); + bmpy::add_dist_params(c, {BMPY_PARAM("probabilities", probabilities), BMPY_PARAM("rates", rates)}); } } diff --git a/nanobind/src/distributions_noncentral.cpp b/nanobind/src/distributions_noncentral.cpp index 0e86d5c363..1aa2d335d4 100644 --- a/nanobind/src/distributions_noncentral.cpp +++ b/nanobind/src/distributions_noncentral.cpp @@ -22,23 +22,27 @@ void bmpy::bind_distributions_noncentral(nb::module_& m) { nb::class_ c(m, "non_central_beta_distribution"); c.def(nb::init(), arg("alpha"), arg("beta"), arg("lambda_")); add_dist_methods(c); + bmpy::add_dist_params(c, {BMPY_PARAM("alpha", alpha), BMPY_PARAM("beta", beta), BMPY_PARAM("lambda_", non_centrality)}); } { using D = bm::non_central_chi_squared_distribution; nb::class_ c(m, "non_central_chi_squared_distribution"); c.def(nb::init(), arg("df"), arg("lambda_")); add_dist_methods(c); + bmpy::add_dist_params(c, {BMPY_PARAM("df", degrees_of_freedom), BMPY_PARAM("lambda_", non_centrality)}); } { using D = bm::non_central_f_distribution; nb::class_ c(m, "non_central_f_distribution"); c.def(nb::init(), arg("df1"), arg("df2"), arg("lambda_")); add_dist_methods(c); + bmpy::add_dist_params(c, {BMPY_PARAM("df1", degrees_of_freedom1), BMPY_PARAM("df2", degrees_of_freedom2), BMPY_PARAM("lambda_", non_centrality)}); } { using D = bm::non_central_t_distribution; nb::class_ c(m, "non_central_t_distribution"); c.def(nb::init(), arg("df"), arg("delta")); add_dist_methods(c); + bmpy::add_dist_params(c, {BMPY_PARAM("df", degrees_of_freedom), BMPY_PARAM("delta", non_centrality)}); } } diff --git a/nanobind/tests/test_smoke.py b/nanobind/tests/test_smoke.py index 5dbb5eb9a9..a1b6d4bf6e 100644 --- a/nanobind/tests/test_smoke.py +++ b/nanobind/tests/test_smoke.py @@ -421,6 +421,7 @@ def test_array_errors_warn_and_give_nan(): ("chi_squared_distribution", (3.5,)), ("exponential_distribution", (2.0,)), ("extreme_value_distribution", ()), + ("generalized_extreme_value_distribution", (0.5, 2.0, 0.1)), ("fisher_f_distribution", (3.0, 7.0)), ("gamma_distribution", (2.5, 1.5)), ("inverse_chi_squared_distribution", (4.0,)), @@ -503,6 +504,23 @@ def test_discrete_distribution(case): assert math.isfinite(dist.mean()) and math.isfinite(dist.variance()) +@pytest.mark.parametrize("case", CONTINUOUS + DISCRETE, ids=lambda c: c[0]) +def test_distribution_repr(case): + # The repr names the constructor's arguments, each a property, and rebuilds the distribution. + name, params = case + dist = getattr(bm, name)(*params) + r = repr(dist) + assert r.startswith(name + "(") + again = eval(r, vars(bm)) + assert repr(again) == r + assert np.allclose(again.cdf(dist.quantile(0.3)), 0.3, rtol=1e-12) or name in {c[0] for c in DISCRETE} + for i, value in enumerate(params): + arg = r[len(name) + 1:-1].split(", ")[i].split("=")[0] if name != "hyperexponential_distribution" else None + if arg is not None and not isinstance(value, (list, tuple)): + got = getattr(dist, arg) + assert (got() if callable(got) else got) == value + + def test_normal_against_mpmath(): dist = bm.normal_distribution(1.0, 2.0) for x in (-3.0, 0.5, 4.0): diff --git a/reporting/performance/generalized_extreme_value_performance.cpp b/reporting/performance/generalized_extreme_value_performance.cpp new file mode 100644 index 0000000000..a6b48a9b7c --- /dev/null +++ b/reporting/performance/generalized_extreme_value_performance.cpp @@ -0,0 +1,109 @@ +// (C) Copyright Nick Thompson 2026. +// Use, modification and distribution are subject to the +// Boost Software License, Version 1.0. (See accompanying file +// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) + +// +// The generalized extreme value distribution's pdf, cdf and quantile at shape xi = 0.2, against the +// Gumbel distribution (extreme_value, which is xi = 0) as the baseline, and its skewness, whose +// moment series is the expensive part of the moments. +// +// Promotion is disabled so that each type is timed in its own arithmetic: by default, +// float is evaluated in double, and double in long double where that is wider. +// + +#include +#include +#include +#include +#include + +using no_promotion = boost::math::policies::policy, boost::math::policies::promote_double>; + +template +std::vector abscissas() +{ + std::mt19937_64 mt {42}; + std::uniform_real_distribution dis(-2, 20); + std::vector v(1024); + for (auto& x : v) + { + x = static_cast(dis(mt)); + } + return v; +} + +template +std::vector probabilities() +{ + std::mt19937_64 mt {42}; + std::uniform_real_distribution dis(0.001, 0.999); + std::vector v(1024); + for (auto& p : v) + { + p = static_cast(dis(mt)); + } + return v; +} + +template +auto gev() +{ + return boost::math::generalized_extreme_value_distribution(Real(0.5), Real(2), Real(0.2)); +} + +template +auto gumbel() +{ + return boost::math::extreme_value_distribution(Real(0.5), Real(2)); +} + +#define BENCH(name, make, fn, points) \ +template \ +void name(benchmark::State& state) \ +{ \ + auto d = make(); \ + auto v = points(); \ + std::size_t i = 0; \ + for (auto _ : state) \ + { \ + benchmark::DoNotOptimize(fn(d, v[i++ & 1023])); \ + } \ +} + +BENCH(gev_pdf, gev, boost::math::pdf, abscissas) +BENCH(gumbel_pdf, gumbel, boost::math::pdf, abscissas) +BENCH(gev_cdf, gev, boost::math::cdf, abscissas) +BENCH(gumbel_cdf, gumbel, boost::math::cdf, abscissas) +BENCH(gev_quantile, gev, boost::math::quantile, probabilities) +BENCH(gumbel_quantile, gumbel, boost::math::quantile, probabilities) + +template +void gev_skewness(benchmark::State& state) +{ + std::mt19937_64 mt {42}; + std::uniform_real_distribution dis(-0.5, 0.3); + std::vector> v; + for (int i = 0; i < 1024; ++i) + { + v.emplace_back(Real(0), Real(1), static_cast(dis(mt))); + } + std::size_t i = 0; + for (auto _ : state) + { + benchmark::DoNotOptimize(boost::math::skewness(v[i++ & 1023])); + } +} + +BENCHMARK_TEMPLATE(gev_pdf, float); +BENCHMARK_TEMPLATE(gev_pdf, double); +BENCHMARK_TEMPLATE(gumbel_pdf, double); +BENCHMARK_TEMPLATE(gev_cdf, float); +BENCHMARK_TEMPLATE(gev_cdf, double); +BENCHMARK_TEMPLATE(gumbel_cdf, double); +BENCHMARK_TEMPLATE(gev_quantile, float); +BENCHMARK_TEMPLATE(gev_quantile, double); +BENCHMARK_TEMPLATE(gumbel_quantile, double); +BENCHMARK_TEMPLATE(gev_skewness, double); + +BENCHMARK_MAIN(); diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 29bcf179aa..f03362f7a9 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -38,12 +38,16 @@ if(HAVE_BOOST_TEST) message(STATUS "Building boost.math with SYCL") + # icpx defaults to -fp-model=fast, which flushes denormals to zero and reassociates arithmetic, + # defeating the error-free transforms and compensated sums the library relies on; test as + # the other compilers do (contraction stays on). + set(CMAKE_CXX_COMPILER "icpx") set(CMAKE_C_COMPILER "icx") enable_testing() - boost_test_jamfile(FILE sycl_jamfile LINK_LIBRARIES Boost::math Boost::assert Boost::concept_check Boost::config Boost::core Boost::integer Boost::lexical_cast Boost::multiprecision Boost::predef Boost::random Boost::throw_exception Boost::unit_test_framework sycl COMPILE_DEFINITIONS BOOST_MATH_ENABLE_SYCL=1 COMPILE_OPTIONS -fsycl ) + boost_test_jamfile(FILE sycl_jamfile LINK_LIBRARIES Boost::math Boost::assert Boost::concept_check Boost::config Boost::core Boost::integer Boost::lexical_cast Boost::multiprecision Boost::predef Boost::random Boost::throw_exception Boost::unit_test_framework sycl COMPILE_DEFINITIONS BOOST_MATH_ENABLE_SYCL=1 COMPILE_OPTIONS -fsycl -fp-model=precise ) elseif (BOOST_MATH_BUILD_MODULE) message(STATUS "Building boost.math as a C++20 module") diff --git a/test/Jamfile.v2 b/test/Jamfile.v2 index c558f5c321..359080af69 100644 --- a/test/Jamfile.v2 +++ b/test/Jamfile.v2 @@ -693,6 +693,7 @@ test-suite distribution_tests : [ run test_dist_overloads.cpp /boost/test//boost_unit_test_framework ] [ run test_exponential_dist.cpp /boost/test//boost_unit_test_framework ] [ run test_extreme_value.cpp /boost/test//boost_unit_test_framework ] + [ run test_generalized_extreme_value.cpp ] [ run test_find_location.cpp pch /boost/test//boost_unit_test_framework ] [ run test_find_scale.cpp pch /boost/test//boost_unit_test_framework ] [ run test_fisher_f.cpp /boost/test//boost_unit_test_framework ] @@ -750,28 +751,28 @@ test-suite distribution_tests : [ run test_logistic_dist.cpp /boost/test//boost_unit_test_framework ] [ run test_logcdf.cpp ] [ run test_lognormal.cpp /boost/test//boost_unit_test_framework ] - [ run test_negative_binomial.cpp /boost/test//boost_unit_test_framework + [ run test_negative_binomial.cpp : # command line : # input files : # requirements TEST_FLOAT intel:off : test_negative_binomial_float ] - [ run test_negative_binomial.cpp /boost/test//boost_unit_test_framework + [ run test_negative_binomial.cpp : # command line : # input files : # requirements TEST_DOUBLE intel:off : test_negative_binomial_double ] - [ run test_negative_binomial.cpp /boost/test//boost_unit_test_framework + [ run test_negative_binomial.cpp : # command line : # input files : # requirements TEST_LDOUBLE intel:off : test_negative_binomial_long_double ] - [ run test_negative_binomial.cpp /boost/test//boost_unit_test_framework + [ run test_negative_binomial.cpp : # command line : # input files : # requirements @@ -1650,6 +1651,7 @@ test-suite concepts : [ run compile_test/test_dist_complement_incl.cpp compile_test_main : : : [ check-target-builds ../config//is_ci_sanitizer_run "Sanitizer CI run" : no ] ] [ run compile_test/test_dist_exponential_incl.cpp compile_test_main : : : [ check-target-builds ../config//is_ci_sanitizer_run "Sanitizer CI run" : no ] ] [ run compile_test/test_dist_extreme_val_incl.cpp compile_test_main : : : [ check-target-builds ../config//is_ci_sanitizer_run "Sanitizer CI run" : no ] ] + [ run compile_test/test_dist_generalized_extreme_value_incl.cpp compile_test_main : : : [ check-target-builds ../config//is_ci_sanitizer_run "Sanitizer CI run" : no ] ] [ run compile_test/test_dist_find_location_incl.cpp compile_test_main : : : [ check-target-builds ../config//is_ci_sanitizer_run "Sanitizer CI run" : no ] ] [ run compile_test/test_dist_find_scale_incl.cpp compile_test_main : : : [ check-target-builds ../config//is_ci_sanitizer_run "Sanitizer CI run" : no ] ] [ run compile_test/test_dist_fisher_f_incl.cpp compile_test_main : : : [ check-target-builds ../config//is_ci_sanitizer_run "Sanitizer CI run" : no ] ] diff --git a/test/compile_test/distribution_concept_check.cpp b/test/compile_test/distribution_concept_check.cpp index 3f5253cd71..0102eaa3fc 100644 --- a/test/compile_test/distribution_concept_check.cpp +++ b/test/compile_test/distribution_concept_check.cpp @@ -27,6 +27,7 @@ void instantiate(RealType) function_requires > >(); function_requires > >(); function_requires > >(); + function_requires > >(); function_requires > >(); function_requires > >(); function_requires > >(); diff --git a/test/compile_test/instantiate.hpp b/test/compile_test/instantiate.hpp index 5508268903..89999e7947 100644 --- a/test/compile_test/instantiate.hpp +++ b/test/compile_test/instantiate.hpp @@ -120,6 +120,7 @@ void instantiate(RealType) function_requires > >(); function_requires > >(); function_requires > >(); + function_requires > >(); function_requires > >(); function_requires > >(); function_requires > >(); @@ -162,6 +163,7 @@ void instantiate(RealType) function_requires > >(); function_requires > >(); function_requires > >(); + function_requires > >(); function_requires > >(); function_requires > >(); function_requires > >(); @@ -199,6 +201,7 @@ void instantiate(RealType) function_requires >(); function_requires >(); function_requires >(); + function_requires >(); function_requires >(); function_requires >(); function_requires >(); diff --git a/test/compile_test/test_dist_generalized_extreme_value_incl.cpp b/test/compile_test/test_dist_generalized_extreme_value_incl.cpp new file mode 100644 index 0000000000..313f88611e --- /dev/null +++ b/test/compile_test/test_dist_generalized_extreme_value_incl.cpp @@ -0,0 +1,25 @@ +// Copyright Nick Thompson 2026. +// Use, modification and distribution are subject to the +// Boost Software License, Version 1.0. (See accompanying file +// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) +// +// Basic sanity check that header +// #includes all the files that it needs to. +// +#include +// +// Note this header includes no other headers, this is +// important if this test is to be meaningful: +// +#include "test_compile_result.hpp" + +void compile_and_link_test() +{ + TEST_DIST_FUNC(generalized_extreme_value) +} + +template class boost::math::generalized_extreme_value_distribution >; +template class boost::math::generalized_extreme_value_distribution >; +#ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS +template class boost::math::generalized_extreme_value_distribution >; +#endif diff --git a/test/math_unit_test.hpp b/test/math_unit_test.hpp index db3776077b..b165a53414 100644 --- a/test/math_unit_test.hpp +++ b/test/math_unit_test.hpp @@ -307,11 +307,11 @@ bool check_absolute_error(PreciseReal expected1, Real computed, Real acceptable_ if (!std::is_integral::value) { BOOST_MATH_ASSERT_MSG(sizeof(PreciseReal) >= sizeof(Real), "The expected number must be computed in higher (or equal) precision than the number being tested."); - BOOST_MATH_ASSERT_MSG(!isnan(expected1), "Expected value cannot be a nan (use CHECK_NAN if this is your intention)."); + BOOST_MATH_ASSERT_MSG(!boost::math::isnan(expected1), "Expected value cannot be a nan (use CHECK_NAN if this is your intention)."); } BOOST_MATH_ASSERT_MSG(acceptable_error > 0, "Error must be > 0."); - if (isnan(computed)) + if (boost::math::isnan(computed)) { std::ios_base::fmtflags f(std::cerr.flags()); std::cerr << std::setprecision(3); @@ -499,8 +499,15 @@ inline int report_errors() #define CHECK_TRUE(X) boost::math::test::check_true((X), __FILE__, __func__, __LINE__) +// Without exceptions (GPU and SYCL builds, -fno-exceptions) errors don't throw, so there is nothing to check; +// as BOOST_MATH_CHECK_THROW does there, the expressions aren't evaluated. +#ifdef BOOST_MATH_NO_EXCEPTIONS +#define CHECK_THROW(x, what) +#define CHECK_NO_THROW(x) x; +#else #define CHECK_THROW(x, what) try{ x; boost::math::test::report_non_throw(__FILE__, __LINE__); }catch(const what&){} catch(...){ boost::math::test::report_incorrect_throw(__FILE__, __LINE__); } #define CHECK_NO_THROW(x) try{ x; }catch(...){ boost::math::test::report_unexpected_throw(__FILE__, __LINE__); } +#endif #endif diff --git a/test/test_dist_deduction_guides.cpp b/test/test_dist_deduction_guides.cpp index b814ee1503..a65d6a1e0e 100644 --- a/test/test_dist_deduction_guides.cpp +++ b/test/test_dist_deduction_guides.cpp @@ -19,6 +19,7 @@ #include #include #include +#include #include //#include //#include @@ -84,6 +85,10 @@ int main() test_deduction_guide(1); test_deduction_guide(1, 1); + test_deduction_guide(0); + test_deduction_guide(0, 1); + test_deduction_guide(0, 1, 0); + test_deduction_guide(1); test_deduction_guide(1); diff --git a/test/test_generalized_extreme_value.cpp b/test/test_generalized_extreme_value.cpp new file mode 100644 index 0000000000..3101ed0fb3 --- /dev/null +++ b/test/test_generalized_extreme_value.cpp @@ -0,0 +1,443 @@ +// Copyright Nick Thompson 2026. +// Distributed under the Boost Software License, Version 1.0. +// (See accompanying file LICENSE_1_0.txt or copy at +// https://www.boost.org/LICENSE_1_0.txt) +// +// References: mpmath at 150 digits, printed to 40, at the double values of the parameters and arguments +// (the generator is in the PR that added this file). Where the tests need more than a few ulps, it is the +// condition number of the function, not the method: +// * the cdf deep in its tails, exp(-t) for large t, magnifies the rounding of (x - mu)/sigma, so the +// tolerances scale with each function's condition number; +// * mu + sigma q cancels where the quantile is near 0 (-0.016 at xi = 1.5, p = 0.25), so the quantile's tolerance +// scales with (|mu| + |x - mu|)/|x|; +// * skewness and kurtosis excess pass through 0 (at xi = -0.28 and -0.15), so they are checked to an +// absolute error, and for |xi| >= 1/4 the kurtosis is a fourth difference of lgamma(1 - k xi). + +#include "math_unit_test.hpp" +#include +#include +#include "test_out_of_range.hpp" +#include +#include +#if __has_include() +#include +#define BOOST_MATH_TEST_GEV_MULTIPRECISION +#endif + +using boost::math::generalized_extreme_value_distribution; + +// mu, sigma, xi, x, pdf, logpdf, cdf, ccdf, logcdf, logccdf +const long double cdf_data[][10] = { + { 0.5, 2, 0, -3, 0.009115765822384945444346020789780975286465L, -4.697749856565675746283731826300868938998L, 0.003168165149053243848719669162880521440672L, 0.9968318348509467561512803308371194785593L, -5.754602676005730436866499704842692370923L, -0.003173194409419889751857546828551054431638L }, + { 0.5, 2, 0, -0.5, 0.1585209605389710880473291448008405244162L, -1.841868451260073456265882909272340139729L, 0.1922956455479649280742866860378975816894L, 0.8077043544520350719257133139621024183106L, -1.648721270700128146848650787814163571654L, -0.2135591853734341440449422035125653857933L }, + { 0.5, 2, 0, 0.25, 0.1824469864020708676088842652042828343706L, -1.701295633626771626246239349269970440641L, 0.3220178007714426731329932585931980462668L, 0.6779821992285573268670067414068019537332L, -1.133148453066826316829007227811793872566L, -0.3886342462115421229163077265709873030708L }, + { 0.5, 2, 0, 1, 0.1787176730860912736876052717082097864081L, -1.721947963631350177662402388436497215372L, 0.4589560693076638184581480769101370726848L, 0.5410439306923361815418519230898629273152L, -0.7788007830714048682451702669783206472968L, -0.614254800673943822670724491688423891635L }, + { 0.5, 2, 0, 4, 0.07302735923472654790400422830811197399479L, -2.61692112401039043609794938012454758409L, 0.8404868737475783667364934360395455344152L, 0.1595131262524216332635065639604544655848L, -0.1737739434504451266807172586663710160147L, -1.835629063889961509423681941292199838984L }, + { 0.5, 2, 0, 30, 1.963931001335069415616541086995321911745e-7L, -15.44314757334629985752113467810391568766L, 0.9999996072137225924461571189333918174662L, 3.927862774075538428810666081825337687368e-7L, -3.927863545481039025566457391195863032792e-7L, -14.75000019639317084567193797879849698995L }, + { 0.5, 2, 1e-20, -3, 0.009115765822384945443841874696649368440024L, -4.697749856565675746339036679777206683281L, 0.003168165149053243848440498834569982501757L, 0.9968318348509467561515595011654300174982L, -5.754602676005730436954617058319030113406L, -0.003173194409419889751577489231521880065215L }, + { 0.5, 2, 1e-20, -0.5, 0.1585209605389710880479932047048435720762L, -1.841868451260073456261693810860715300143L, 0.192295645547964928073890383636550153991L, 0.807704354452035071926109616363449846009L, -1.648721270700128146850711689402538731724L, -0.2135591853734341440444515507147154989699L }, + { 0.5, 2, 1e-20, 0.25, 0.1824469864020708676091104260830661932527L, -1.701295633626771626244999751492866286515L, 0.3220178007714426731329647512515727226953L, 0.6779821992285573268670352487484272773047L, -1.133148453066826316829095755034689718367L, -0.3886342462115421229162656792404559695394L }, + { 0.5, 2, 1e-20, 1, 0.1787176730860912736871708313409126657638L, -1.72194796363135017766483326368120702926L, 0.4589560693076638184580363783644582656449L, 0.5410439306923361815419636216355417343551L, -0.7788007830714048682454136422230304613048L, -0.6142548006739438226705180416400223243136L }, + { 0.5, 2, 1e-20, 4, 0.07302735923472654790365016147306552690034L, -2.616921124010390436102797793633632524827L, 0.8404868737475783667342569731629820340084L, 0.1595131262524216332657430268370179659916L, -0.1737739434504451266833781721754559568711L, -1.83562906388996150940966138430666117859L }, + { 0.5, 2, 1e-20, 30, 1.963931001335069417463249117543877889429e-7L, -15.44314757334629985658082260538182204619L, 0.9999996072137225924461566916556533395052L, 3.927862774075538433083443466604947852275e-7L, -3.927863545481039029839236454264455565347e-7L, -14.75000019639317084458412569243743621741L }, + { 0.5, 2, -1e-9, -3, 0.009115765872799554771946080968262620026985L, -4.697749851035190412359091169626025457627L, 0.003168165176970276748828468916096131826271L, 0.9968318348230297232511715310839038681737L, -5.754602667193995106259567173008525775103L, -0.003173194437425649524441523294602566123376L }, + { 0.5, 2, -1e-9, -0.5, 0.1585209604725650976744218719492888455378L, -1.841868451678983297343277604902220403243L, 0.1922956455875951681974261917104534818515L, 0.8077043544124048318025738082895465181485L, -1.648721270494037988092712111127424194838L, -0.2135591854224994238167688548412075966868L }, + { 0.5, 2, -1e-9, 0.25, 0.182446986379454779733353359398335601163L, -1.701295633750731403948963542528035519475L, 0.3220178007742934072952895375955548586645L, 0.6779821992257065927047104624044451413355L, -1.133148453057973594540194954464039676811L, -0.3886342462157468559691017126026166415997L }, + { 0.5, 2, -1e-9, 1, 0.178717673129535310427908761610173438766L, -1.721947963388262653161703324233783748909L, 0.4589560693188336730278518189664221151632L, 0.5410439306811663269721481810335778848368L, -0.7788007830470673437705128830706037529678L, -0.6142548006945888275144638144559958666686L }, + { 0.5, 2, -1e-9, 4, 0.07302735927013323140638808443280455603325L, -2.616921123525549085337949864699606090224L, 0.8404868739712246545122901543310216909817L, 0.1595131260287753454877098456689783090183L, -0.17377394318435377566550942297392605831L, -1.835629065292017209709248439735100776865L }, + { 0.5, 2, -1e-9, 30, 1.963930816664273156023001573027157430863e-7L, -15.4431476673775080906334661135206673916L, 0.9999996072137653202181010821144661558863L, 3.927862346797818989178855338441136576943e-7L, -3.927863118203151757053016200034791262688e-7L, -14.7500003051744005514612834025113567628L }, + { 0.5, 2, 1e-9, -3, 0.009115765771970336145624700217464420508732L, -4.697749862096161107626641378518741212887L, 0.003168165121136211086720674613104760872243L, 0.9968318348788637889132793253868952391278L, -5.75460268481746580152711779888655526371L, -0.003173194381414130118606643465339249650594L }, + { 0.5, 2, 1e-9, -0.5, 0.1585209606053770784750314195231078862963L, -1.84186845084116361501830958938882670287L, 0.1922956455083346879279406427428083575554L, 0.8077043544916653120720593572571916424446L, -1.648721270906218305767744173580603150022L, -0.213559185324368864246791478672148497165L }, + { 0.5, 2, 1e-9, 0.25, 0.1824469864246869554900251411162551697432L, -1.701295633502811848528132688286390169526L, 0.3220178007685919389702248226276367206498L, 0.6779821992314080610297751773723632793502L, -1.133148453075679039119364116767366211481L, -0.388634246207337389862835005334260806701L }, + { 0.5, 2, 1e-9, 1, 0.1787176730426472369684846225640671376264L, -1.721947963874437702103666158788601709202L, 0.4589560692964939638920904083634490935538L, 0.5410439307035060361079095916365509064462L, -0.7788007830957423927124756903687619088155L, -0.6142548006532988178341503453350002029345L }, + { 0.5, 2, 1e-9, 4, 0.07302735919931986439709865733298120907751L, -2.616921124495231787154938124508301140614L, 0.8404868735239320791995900204472480338069L, 0.1595131264760679208004099795527519661931L, -0.1737739437165364774824976566509613042549L, -1.835629062487905812601516027788899559118L }, + { 0.5, 2, 1e-9, 30, 1.963931186005879265734471889093366625678e-7L, -15.44314747931509354621469424527979978518L, 0.999999607213679864670405489814773606136L, 3.927863201353295945101852263938640014138e-7L, -3.927863972758964370772486657174537095069e-7L, -14.75000008761194327924907972277214799304L }, + { 0.5, 2, 0.05, -3, 0.006654207329959662526934006477968915074622L, -5.012505943175803363223445387272214509505L, 0.001945434707953234552092815540172502183168L, 0.9980545652920467654479071844598274978168L, -6.242269826651158210658094099245287635833L, -0.00194732952394737760557329759096809338509L }, + { 0.5, 2, 0.05, -0.5, 0.1619111752996693505827179422710667011729L, -1.820707394740831011982784453049470024416L, 0.1902846477536663758749416549401691954994L, 0.8097153522463336241250583450598308045006L, -1.659234181850973087825942269518598513973L, -0.2110725100569042567305917116110946846935L }, + { 0.5, 2, 0.05, 0.25, 0.1835848470423718787177487583081779818812L, -1.695078336650691841412114949055366012756L, 0.3218746711043940370184449592529980726989L, 0.6781253288956059629815550407470019273011L, -1.1335930293762498514571924436629684101L, -0.3884231572460850972272443880540950963379L }, + { 0.5, 2, 0.05, 1, 0.1765716300692055618376724509591067746702L, -1.734028648851825495277856801880003561455L, 0.4584020927196983646983517538704278540592L, 0.5415979072803016353016482461295721459408L, -0.7800085483221799657234567565864361429962L, -0.6132314213496210496449898360512151050468L }, + { 0.5, 2, 0.05, 4, 0.07125611704849486727158939803010069722299L, -2.641474610206353167989633676116143473334L, 0.8295961473591191058350015746177976544793L, 0.1704038526408808941649984253822023455207L, -0.1868162660516636254194341504841769437842L, -1.769584055447669281973206403773798709673L }, + { 0.5, 2, 0.05, 30, 4.576590121622874632575130895371787826148e-6L, -12.29455635205656301223688087145681232746L, 0.9999840962228616264473566472665571532576L, 1.590377713837355264335273344284674240573e-5L, -1.590390360477805028881394339502371816681e-5L, -11.04895392107746622240188923278841501568L }, + { 0.5, 2, -0.05, -3, 0.0116522972401220260452493356780315964923L, -4.452251930416117493176344205665701826687L, 0.004734624088379052327935505548823451317078L, 0.9952653759116209476720644944511765486829L, -5.352852945489512195361591654164601106664L, -0.004745867925273727842992985544106665896138L }, + { 0.5, 2, -0.05, -0.5, 0.1552677890261829157254418743604834157109L, -1.86260398163328395446559383287104030657L, 0.1942486110094597440911993454966736347911L, 0.8057513889905402559088006545033263652089L, -1.638616440290397162630580347606087362705L, -0.2159800344400579524535708917838999447462L }, + { 0.5, 2, -0.05, 0.25, 0.181323151227319643223971569411501005637L, -1.707474473690776936044334692444143202217L, 0.3221597500187632053101022398452834952805L, 0.6778402499812367946898977601547165047195L, -1.132707738392917036539335266553895759585L, -0.3888436383010770192104506812337850194648L }, + { 0.5, 2, -0.05, 1, 0.1809166811361984369704735047807703523638L, -1.709718678912861502344432957419012997998L, 0.4595191619005060828100913912866177786501L, 0.5404808380994939171899086087133822213499L, -0.7775746364225739497035679822722962193043L, -0.6152960947061666546825125553326396886438L }, + { 0.5, 2, -0.05, 4, 0.07478619311646186860222396651672830589823L, -2.593121995043949955032873487064120053411L, 0.8519749640523185531180141938235932211261L, 0.1480250359476814468819858061764067788739L, -0.1601981374996854657648926641497742859958L, -1.910373857715608960778425274430996387585L }, + { 0.5, 2, -0.05, 30, 4.596495628665318784510202472933638098048e-12L, -26.10572692262107364754515967239754052211L, 0.9999999999975868397949477963434007444284L, 2.413160205052203656599255571631276795675e-12L, -2.413160205055115327686884052649638212485e-12L, -26.7500839390103805297543906676252635749L }, + { 0.5, 2, 0.125, -3, 0.003423246589701594419900482066542693821116L, -5.677165882969006728455454653911608516139L, 0.0007422982779016125799840396566571464186797L, 0.9992577017220983874200159603433428535813L, -7.2057594037927936L, -0.0007425739176813640686958440054640087028603L }, + { 0.5, 2, 0.125, -0.5, 0.1672727062307847691187462917250953461827L, -1.788129826971705694097237563124980650824L, 0.1871529023412981934236701852128733036644L, 0.8128470976587018065763298147871266963356L, -1.675829336649900929736320682822740435909L, -0.2072122588825521358688092544051501334343L }, + { 0.5, 2, 0.125, 0.25, 0.1853185158707079340417764264213095687752L, -1.685679226933769823085801112302159614852L, 0.3216577267813451408584717373817422758756L, 0.6783422732186548591415282626182577241244L, -1.134267259087077031136514593991437472465L, -0.3881032906682759843935688114900746473114L }, + { 0.5, 2, 0.125, 1, 0.1734477527072901158512380207091809130925L, -1.751878862079889173325343037941231293253L, 0.4575876840435851385758078160443173339454L, 0.5424123159564148614241921839556826660546L, -0.7817867535191606685688570481121052483518L, -0.6117288362044184134725997865669490731972L }, + { 0.5, 2, 0.125, 4, 0.06863052259648018182869424109167513063265L, -2.679017907396134254294770921115472364825L, 0.8142897347074998526170434672704840171996L, 0.1857102652925001473829565327295159828004L, -0.2054390368669100216143901588846425750722L, -1.683567533247151993873221105493973708099L }, + { 0.5, 2, 0.125, 30, 4.110037028923576548359531029915430181689e-5L, -10.09949342703221363549386646401175817482L, 0.9997662143139766060126797142254146386006L, 0.0002337856860233939873202857745853613993746L, -0.0002338130181568794225181430950421601604372L, -8.361105733968211082512536218011388359181L }, + { 0.5, 2, 0.13, -3, 0.003238608780817604308032961692976858258538L, -5.732611430028139465850306873793424785492L, 0.0006870172638822821523678523743578983473903L, 0.9993129827361177178476321476256421016526L, -7.283151136680250121578402982594014422121L, -0.0006872533683874924330563580122977629701224L }, + { 0.5, 2, 0.13, -0.5, 0.167642582979750668189820404169765257474L, -1.785921048147577612224991852149429323638L, 0.186938984860475730721935780976793185975L, 0.813061015139524269278064219023206814025L, -1.676972999538390460021589330236193440741L, -0.2069491228782811202391501773295148965376L }, + { 0.5, 2, 0.13, 0.25, 0.1854352593863591939810156732461367276423L, -1.685049463854870883476988933156621950072L, 0.3216431666550400036284121905812256791238L, 0.6783568333449599963715878094187743208762L, -1.134312526011019406511383155627553356037L, -0.3880818266218485114822986459015728619994L }, + { 0.5, 2, 0.13, 1, 0.173243430462994115239555769316327097926L, -1.753057561036600643732551237372863172629L, 0.4575340801157133252464811627294380131082L, 0.5424659198842866747535188372705619868918L, -0.7819039049847524964377866345172480098276L, -0.611630016030606195189945313598941525237L }, + { 0.5, 2, 0.13, 4, 0.06845798596980971827317031124342798230274L, -2.681535065392460035884664670868056662852L, 0.813310520058397836922121976181850860999L, 0.186689479941602163077878023818149139001L, -0.2066422988649100421747385705977984967315L, -1.678308577415702260218934747197293895355L }, + { 0.5, 2, 0.13, 30, 4.537717623727485757201078903810861975631e-5L, -10.00050130537729555640686691536980221853L, 0.9997351891111566976157126505655734167769L, 0.0002648108888433023842873494344265832231021L, -0.0002648459574378944841479534482292109445275L, -8.236494613737416399497571797131769469261L }, + { 0.5, 2, -0.3, -3, 0.02257883218224408674428728479627418398383L, -3.79074244071926629082213152250139356226L, 0.01686948545248619375262787356543226767867L, 0.9831305145475138062473721264345677323213L, -4.082248883631195877498228653947753613791L, -0.01701339598251505476486047576402773046538L }, + { 0.5, 2, -0.3, -0.5, 0.1407963398662055892525447031129239136436L, -1.960440830864504084642685787505986146457L, 0.2032324564965004859201866072448043529658L, 0.7967675435034995140798133927551956470342L, -1.59340484917992907506999272775344169187L, -0.227192307096227335984562664303698179324L }, + { 0.5, 2, -0.3, 0.25, 0.1759055476293873263072756029607300024377L, -1.737808089194809348779680153454023765511L, 0.3228523505680641335944189777181274753936L, 0.6771476494319358664055810222818725246064L, -1.130560179254535433596138358277503983987L, -0.3898659361532131300854343785533963854326L }, + { 0.5, 2, -0.3, 1, 0.1927795672983714006835847556170475107805L, -1.646207881164722508291043440915108869867L, 0.4624807071474759097728591303018553747857L, 0.5375192928525240902271408696981446252143L, -0.7711504371754495264611890107350469005594L, -0.6207906259453672883184693821818092690651L }, + { 0.5, 2, -0.3, 4, 0.08096038522226132802083470239305868591372L, -2.513795315274196404958072065994092601609L, 0.9197803753929363396011351962134600433416L, 0.08021962460706366039886480378653995665844L, -0.08362035983676079918029272072320952064178L, -2.522987098194140438626166073892907539789L }, + { 0.5, 2, 0.3, -3, 8.059704992362896995119811007517873715506e-5L, -9.426048510565349992020177118882759848275L, 6.402576600575467769883962668923564384817e-6L, 0.9999935974233994245322301160373310764356L, -11.95881005477788658284444167501990557427L, -6.402597097156518221636218314392215847589e-6L }, + { 0.5, 2, 0.3, -0.5, 0.1812503225698147173179925232617782002554L, -1.707876205553407075855624216073082303607L, 0.1792492973447090255388108136039045090739L, 0.8207507026552909744611891863960954909261L, -1.71897771948381971532460359928267808018L, -0.1975358664871359474170686040679475679825L }, + { 0.5, 2, 0.3, 0.25, 0.1894940236404238306364162765021104566183L, -1.663397792474526321406233515749467943149L, 0.3211407214070816380060702435413827912124L, 0.6788592785929183619939297564586172087876L, -1.135875867468771317232919576205998713349L, -0.3873414209244915061172579286398172204468L }, + { 0.5, 2, 0.3, 1, 0.1665725681295165996385848051549330081597L, -1.79232421989259561589974139903168544228L, 0.4557606462631520925261687698121302605051L, 0.5442393537368479074738312301878697394949L, -0.7857875058209371193944533846408445480066L, -0.6083661404058928633051700276338946733242L }, + { 0.5, 2, 0.3, 4, 0.06286582155301994052150842939996142150539L, -2.766752640494270141607384901313505234844L, 0.7827334860527895458826706624083477548845L, 0.2172665139472104541173293375916522451155L, -0.2449630163437000697794851493209835780066L, -1.526630504159534364420939159431703465001L }, + { 0.5, 2, 0.3, 30, 0.0003273580834484173214067780354743419595302L, -8.024455929786768161735886630208929083997L, 0.996441826975464168416324204254454132273L, 0.003558173024535831583675795745545867727038L, -0.003564518378555365734021672304781503586974L, -5.638508061203108378546491617990135356678L }, + { 0.5, 2, -0.7, -3, 0.03065313768572185299711007995410933359305L, -3.485020250364091045199580492854250520331L, 0.04351602431599833910256362151992301784644L, 0.9564839756840016608974363784800769821536L, -3.134626033640518737066247611835648581181L, -0.04449124331691377264478415694517319717866L }, + { 0.5, 2, -0.7, -0.5, 0.1224780659178903405081370654853399188579L, -2.099823318760895943245804892832088841208L, 0.2153927028350637145526002290527971174053L, 0.7846072971649362854473997709472028825947L, -1.535292392108238402871228881400685712766L, -0.2425719447593194681492008503702743807618L }, + { 0.5, 2, -0.7, 0.25, 0.167880318057502716578366363708646221341L, -1.784503945965868755929679208375181402125L, 0.3239046900990341450574093655119228912695L, 0.6760953099009658549425906344880771087305L, -1.127305972826224354630816284263039890692L, -0.3914212218997065966819085426931104755151L }, + { 0.5, 2, -0.7, 1, 0.2153903268677851220270401478970179674622L, -1.535303423031532923288919312610236116888L, 0.4678015340166633970763864726217391180225L, 0.5321984659833366029236135273782608819775L, -0.7597111456226778590349280458738389255287L, -0.6307388028490597649920765035511791995061L }, + { 0.5, 2, 0.9, -0.5, 0.2530658734606315277479524224191294277061L, -1.374105454735238240559972359301174861478L, 0.1432650870655338913562264826490839940841L, 0.8567349129344661086437735173509160059159L, -1.943058609103825017156942536489297314252L, -0.1546267280300151672391999106680229972773L }, + { 0.5, 2, 0.9, 0.25, 0.2053591919323515073717253479034976085078L, -1.582994677310361952163458363526016508605L, 0.3192433025902280900330646960560710881853L, 0.6807566974097719099669353039439289118147L, -1.141801762863611938061786002862994555557L, -0.3845503091956913171670906007042623815087L }, + { 0.5, 2, 0.9, 1, 0.1466505062859526950907301644183770990309L, -1.919703031222647272046695408322428608797L, 0.450172165222532504053136402084222020012L, 0.549827834777467495946863597915777979988L, -0.7981251800030224208698105134477822631855L, -0.59815007743634220782812536170051942518L }, + { 0.5, 2, 0.9, 4, 0.04785650470556758324311619753289467322999L, -3.039548230811817410079180781740187256459L, 0.7049646710996157775249081280210048415821L, 0.2950353289003842224750918719789951584179L, -0.3496075893409509957988211565573886561168L, -1.220660170828498612057196401644921029177L }, + { 0.5, 2, 0.9, 30, 0.00173336156054610425360097856463650573284L, -6.357692657265018351117128996627085044189L, 0.9491997737464478260487518121425344402256L, 0.05080022625355217395124818785746555977442L, -0.05213599278291702573343607369694041681873L, -2.979854470597536862470282274060356423935L }, + { 0.5, 2, 2, 0.25, 0.2426040447120236162512332757790964753174L, -1.416324610261525447276701173471350332115L, 0.3151518986722023896141958639310454513001L, 0.6848481013277976103858041360689545486999L, -1.154700538379251529018297561003914911295L, -0.3785582151964799813383119465458937224332L }, + { 0.5, 2, 2, 1, 0.1202910057881380141120158758405166072322L, -2.117841423649917915116679819556664070255L, 0.4419773772406739579865398211578859656963L, 0.5580226227593260420134601788421140343037L, -0.816496580927726032732428024901963797322L, -0.5833557748432169237414957712355888063703L }, + { 0.5, 2, 2, 4, 0.03269059698142522532567882310927362868032L, -3.420667796515388102411015891441721522761L, 0.6241250557782609608783975539903324977094L, 0.3758749442217390391216024460096675022906L, -0.4714045207910316829338962414032326928566L, -0.9784987860692273495936509450915494075046L }, + { 0.5, 2, 2, 30, 0.002476748424901686938892812180090879881433L, -6.000808698065031280290777657414978377417L, 0.8343757053045658053987031428930575316221L, 0.1656242946954341946012968571069424683779L, -0.1810714920850370618723100630056445419812L, -1.798033340693379056209812357545670028983L }, + { 0.5, 2, -2, -3, 0.02825439601041371508207652025113584539374L, -3.56650622250772491930639438396616110654L, 0.1198732501037620330747323068407013287011L, 0.8801267498962379669252676931592986712989L, -2.121320343559642573202533086314547117855L, -0.1276893479087839697383530966932791042295L }, + { 0.5, 2, -2, -0.5, 0.08595474576918094590790089441973243385282L, -2.453934333213013012927536906396962930683L, 0.2431167344342142108048623204999460644609L, 0.7568832655657857891951376795000539355391L, -1.41421356237309504880168872420969807857L, -0.278546244085973154673030816903783655471L }, + { 0.5, 2, -2, 0.25, 0.1462039162679206421758576956058319087134L, -1.922752944966945035504966500978731937483L, 0.3269218953517579112720873007606364511319L, 0.6730781046482420887279126992393635488681L, -1.11803398874989484820458683436563811772L, -0.3958939016161662776988365565973799879904L }, + { 0.5, 2, -2, 1, 0.3486522152763511488462237401551361482114L, -1.053680371466520179109460422833937323323L, 0.4930686913952397878459573806151774752358L, 0.5069313086047602121540426193848225247642L, -0.7071067811865475244008443621048490392848L, -0.6793797705566755265109604011331633483493L }, +}; + +// mu, sigma, xi, p, quantile(p), quantile(complement(p)) +const long double quantile_data[][6] = { + { 0.5, 2, 0, 1e-300, -12.5756298398083137183963288862606664562L, 1382.051055796427410360676689140201005817L }, + { 0.5, 2, 0, 1e-10, -5.773235076484002964477948553537579912083L, 46.55170185978091360749126779238166592248L }, + { 0.5, 2, 0, 0.001, -3.365289467832130976369869903602085211958L, 14.31451014104743295818180028793985985517L }, + { 0.5, 2, 0, 0.25, -0.1532685199565619648095859264510141972425L, 2.991798647414476396761561605799278127912L }, + { 0.5, 2, 0, 0.5, 1.233025841163328654024878316465338938909L, 1.233025841163328654024878316465338938909L }, + { 0.5, 2, 0, 0.75, 2.991798647414476396761561605799278127912L, -0.1532685199565619648095859264510141972425L }, + { 0.5, 2, 0, 0.999, 14.31451014104743122259023912057020651812L, -3.365289467832130725242925488683664408872L }, + { 0.5, 2, 0, 0.99999999, 37.34136146785521232696297750383352136488L, -5.326947973310027150864118532956906440812L }, + { 0.5, 2, 1e-9, 1e-300, -12.57562979706528988462358895414425749839L, 1382.05153296736722607948550779053034854L }, + { 0.5, 2, 1e-9, 1e-10, -5.77323506664563339355703516054520922164L, 46.55170238997072872238572308430781231717L }, + { 0.5, 2, 1e-9, 0.001, -3.365289464097015311242587306799330065802L, 14.31451018875760567730617218391453783919L }, + { 0.5, 2, 1e-9, 0.25, -0.1532685198498720250296428667953639746252L, 2.991798648966741522220370573709271182235L }, + { 0.5, 2, 1e-9, 0.5, 1.23302584129766037499459114420234923096L, 1.23302584129766037499459114420234923096L }, + { 0.5, 2, 1e-9, 0.75, 2.991798648966741522220370573709271182235L, -0.1532685198498720250296428667953639746252L }, + { 0.5, 2, 1e-9, 0.999, 14.31451018875760394171459902837123186796L, -3.365289464097015060115643377220075461646L }, + { 0.5, 2, 1e-9, 0.99999999, 37.34136180717669311176557459591262533638L, -5.326947964821696488192169976806530503024L }, + { 0.5, 2, -1e-9, 1e-300, -12.57562988255133773846638851920848108813L, 1382.050578625707339929835571189279739489L }, + { 0.5, 2, -1e-9, 1e-10, -5.773235086322372555971663663944269955237L, 46.55170132959110663131118162134251665025L }, + { 0.5, 2, -1e-9, 0.001, -3.36528947156724664630958691735551269533L, 14.31451009333726045875498283928558450044L }, + { 0.5, 2, -1e-9, 0.25, -0.1532685200632519046127613791272932954428L, 2.991798645862211272592063350702808937606L }, + { 0.5, 2, -1e-9, 0.5, 1.233025841028996933087988362977232413154L, 1.233025841028996933087988362977232413154L }, + { 0.5, 2, -1e-9, 0.75, 2.991798645862211272592063350702808937606L, -0.1532685200632519046127613791272932954428L }, + { 0.5, 2, -1e-9, 0.999, 14.31451009333725872316343366008950099219L, -3.365289471567246395182642017097924755325L }, + { 0.5, 2, -1e-9, 0.99999999, 37.34136112853373570918213062556195007658L, -5.326947981798357830023087488999341876338L }, + { 0.5, 2, 0.1, 1e-300, -9.098492146372804608864205359432952262754L, 2.00000000000000755311483837239514148787e+31L }, + { 0.5, 2, 0.1, 1e-10, -4.884676695574487760921707924416759610688L, 180.4999999990000148431395292428273467397L }, + { 0.5, 2, 0.1, 0.001, -3.014721011455450203818984787388703685023L, 20.40325025512648903682365669916080180916L }, + { 0.5, 2, 0.1, 0.25, -0.1427147455436507519235484201451701256752L, 3.153677616531757077682086383344611767531L }, + { 0.5, 2, 0.1, 0.5, 1.246624642471412121321222799267411548132L, 1.246624642471412121321222799267411548132L }, + { 0.5, 2, 0.1, 0.75, 3.153677616531757077682086383344611767531L, -0.1427147455436507519235484201451701256752L }, + { 0.5, 2, 0.1, 0.999, 20.40325025512648557403643640179963007554L, -3.014721011455449996824097776356218705496L }, + { 0.5, 2, 0.1, 0.99999999, 106.6914687695347467102837941353157503273L, -4.554879256576813897134283068908618375508L }, + { 0.5, 2, -0.1, 1e-300, -17.95596288816075050927756797340028646439L, 20.49999999999999888977697537482352120613L }, + { 0.5, 2, -0.1, 1e-10, -6.868535862554625238918763399218445721122L, 18.49999999998999924915069534853500923024L }, + { 0.5, 2, -0.1, 0.001, -3.764072223342769747986519854910516424594L, 10.47575391872718881130071752413829357759L }, + { 0.5, 2, -0.1, 0.25, -0.164054630693307353519068190091098963657L, 2.842822795884769008582543984252495089731L }, + { 0.5, 2, -0.1, 0.5, 1.219755290644205266101074759104846547934L, 1.219755290644205266101074759104846547934L }, + { 0.5, 2, -0.1, 0.75, 2.842822795884769008582543984252495089731L, -0.164054630693307353519068190091098963657L }, + { 0.5, 2, -0.1, 0.999, 10.47575391872718794140087223803001241633L, -3.764072223342769443318404029362227341775L }, + { 0.5, 2, -0.1, 0.99999999, 17.330213611900137860993975842357316173L, -6.264588046304391117894908529045757023257L }, + { 0.5, 2, 0.5, 1e-300, -3.347808106759349922301748343690987107433L, 3.999999999999999949881816329582481570545e+150L }, + { 0.5, 2, 0.5, 1e-10, -2.666410670026779358283931425752205794977L, 3.999964999899999927132688700516086050025e+5L }, + { 0.5, 2, 0.5, 0.001, -1.978081067593499220751919467145709898194L, 122.9594744018793442987725313740319406105L }, + { 0.5, 2, 0.5, 0.25, -0.1027127988479238291139886635884415211394L, 3.957677382973556235894851336295208920213L }, + { 0.5, 2, 0.5, 0.5, 1.304489635145799179431213144380886890267L, 1.304489635145799179431213144380886890267L }, + { 0.5, 2, 0.5, 0.75, 3.957677382973556235894851336295208920213L, -0.1027127988479238291139886635884415211394L }, + { 0.5, 2, 0.5, 0.999, 122.9594744018792894282733809833461204596L, -1.978081067593499125203206681530880286184L }, + { 0.5, 2, 0.5, 0.99999999, 39996.49979950481432923013743236994190557L, -2.568018796303575793907240515617076686281L }, + { 0.5, 2, -0.5, 1e-300, -100.6304353951386395706955318473740513189L, 4.5L }, + { 0.5, 2, -0.5, 1e-10, -14.69410364875232481508473000078024131598L, 4.499959999999998999999271301886968726567L }, + { 0.5, 2, -0.5, 0.001, -6.013043539513863941419591712727294057535L, 4.373477253676121397844184180029057779858L }, + { 0.5, 2, -0.5, 0.25, -0.2096400900618987640462773058387985509895L, 2.354559914789393416226967493531125213178L }, + { 0.5, 2, -0.5, 0.5, 1.169781555369208974587341420419195809478L, 1.169781555369208974587341420419195809478L }, + { 0.5, 2, -0.5, 0.75, 2.354559914789393416226967493531125213178L, -0.2096400900618987640462773058387985509895L }, + { 0.5, 2, -0.5, 0.999, 4.373477253676121342946231476168160149408L, -6.013043539513863281392466567948484653939L }, + { 0.5, 2, -0.5, 0.99999999, 4.499599999997995048133242721454662269415L, -12.66772820797328732427871819361313767496L }, + { 0.5, 2, 1.5, 1e-200, -0.8331984152353507655021985957119811372087L, 1.333333333333333369132808535167814942226e+300L }, + { 0.5, 2, 1.5, 1e-10, -0.8212658917939733770375279133898950683426L, 1.333333333233332427135188704159432109032e+15L }, + { 0.5, 2, 1.5, 0.001, -0.7598932228613076226472790168569414439187L, 42131.24470703879474650865711481124332923L }, + { 0.5, 2, 1.5, 0.25, -0.01645843373705890274069126157240972584411L, 7.807776763674368551664252384128167313212L }, + { 0.5, 2, 1.5, 0.5, 1.477137790209089876562044068057611286388L, 1.477137790209089876562044068057611286388L }, + { 0.5, 2, 1.5, 0.75, 7.807776763674368551664252384128167313212L, -0.01645843373705890274069126157240972584411L }, + { 0.5, 2, 1.5, 0.999, 42131.24470703873990344933112866136953218L, -0.759893222861307608815186111606538485309L }, + { 0.5, 2, 1.5, 0.99999999, 1.333333313282981533418812268220497771132e+12L, -0.81646857107861184640351984211780469705L }, + { 0.5, 2, -1.5, 1e-300, -24205.34400274864829699997496298744631796L, 1.833333333333333333333333333333333333333L }, + { 0.5, 2, -1.5, 1e-10, -145.4868564499990778008788306250011297411L, 1.833333333333331999999999899999927127689L }, + { 0.5, 2, -1.5, 0.001, -22.37384400274864818889387098884254541468L, 1.83329113798336568488968167869093474847L }, + { 0.5, 2, -1.5, 0.25, -0.3429824999189946874784489084545169083132L, 1.627598449988495743380535845506960621981L }, + { 0.5, 2, -1.5, 0.5, 1.06388949148514695405338186368306265767L, 1.06388949148514695405338186368306265767L }, + { 0.5, 2, -1.5, 0.75, 1.627598449988495743380535845506960621981L, -0.3429824999189946874784489084545169083132L }, + { 0.5, 2, -1.5, 0.999, 1.833291137983365684834756258697668850538L, -22.3738440027486436295880129985934156298L }, + { 0.5, 2, -1.5, 0.99999999, 1.833333333331999999979950481231931414678L, -103.5804134439192277068503656272328464734L }, +}; + +// xi, mean, variance, skewness, kurtosis_excess, mode, median, entropy, for mu = 0.5, sigma = 2 +const long double moment_data[][8] = { + { -3, -2.833333333333333333333333333333333333333L, 304.0L, -19.58485896516058864744704924603960276765L, 1002.531855955678670360110803324099722992L, 1.166666666666666666666666666666666666667L, 0.9446502320073803468540976115921796372292L, 0.5387158507568795882042079412933717059912L }, + { -1, 0.5L, 4.0L, -2.0L, 6.0L, 2.5L, 1.113705638880109381165535757083646863849L, 1.693147180559945309417232121458176568076L }, + { -0.7, 0.7610321918477416120986094092934992813022L, 3.400367327861812827006952124923533645419L, -1.160392728202874330200516009502937526774L, 1.698453497221796334383449480790908683616L, 2.127109656428759269757252685159577858255L, 1.146544828562143746224956954046140776713L, 1.866311880030405193232710602402922212832L }, + { -0.3, 1.183528691291485424098573010753824634699L, 3.913853269298499701383800112689952627786L, -0.06874209942096706080288294973113119974362L, -0.289399162107949706627675781165848166877L, 1.176510388062401562284821150299102160269L, 1.194163629688140530678864252817002609567L, 2.097198145991018318250171797995864498666L }, + { -0.13, 1.424131542591410291946229744699725656592L, 4.952896477265949926405670179598615678555L, 0.5136871089846431160319333561076186998656L, 0.2858096471132615891399975041498497692222L, 0.7760180707317563468194075597259850477742L, 1.215836798220134140775267400701285930257L, 2.195324809024278895581545154437864191538L }, + { -0.12, 1.440163573101605578076943181400254978206L, 5.044061098077296378047764932023288965371L, 0.5540021135248004886617811879635306797331L, 0.3716301414765079931716820994962564984678L, 0.7537157672471984893263687435540211023558L, 1.217139783552452963126246626368292989835L, 2.201096965673294229314315246122693198937L }, + { -0.01, 1.634829761698792572529680260212485936139L, 6.412196648010445039890101186741231605832L, 1.081073759807266273858956528275016526025L, 2.125445886586571376929642583102013399481L, 0.5200996616483335806178703717697134542996L, 1.231684163594858762816687075275816643953L, 2.264590688812462841297521942887004858016L }, + { -1e-6, 1.654429351692890021456710465409222852749L, 6.57971905513086479986738955938859978589L, 1.139541132804515740836640856317379221379L, 2.399971082464274566936674468358194583232L, 0.5000020000009999996665756627996313747545L, 1.2330257068316241121591228516200509969L, 2.27036226824581326849090972512720938759L }, + { -1e-12, 1.654431329801087609222370050051952694087L, 6.579736267375693436975510933807444342468L, 1.139547099398682045080138208052414725792L, 2.399999999971082311659352991890256472674L, 0.5000000000020000000000009999597732949251L, 1.233025841163194322303925031497233263303L, 2.27036284546090095435884267869158222909L }, + { 1e-12, 1.654431329805043833203681940193961353626L, 6.579736267410118054803904146074174922329L, 1.139547099410615269905472390300595465804L, 2.400000000028917688340952240175754611781L, 0.4999999999980000000000010000402267050748L, 1.233025841163462985745831634256318863418L, 2.270362845462055385688645744389575769146L }, + { 1e-6, 1.654433307916871337273665367728207595283L, 6.579753479748693405330144223354347365232L, 1.139553066029341147467124756032196087804L, 2.400028917840957499080165724640479681121L, 0.4999980000010000003334233370193610725493L, 1.233025975495066018764882685290321447262L, 2.270363422677143071556578697953948610645L }, + { 0.01, 1.674395928821558387244769681014319615333L, 6.756655174188248780404314727298438625328L, 1.200478519762556269331078472919269293561L, 2.705134978371804552738435844656986437332L, 0.48010032835166642105707054463640603916L, 1.234370801021428012906777726473035063745L, 2.276135002110493498749966480194153140219L }, + { 0.12, 1.921796457795325141504033129127197829462L, 9.567328041748161760308999329794170737815L, 2.129206380223521923200713361070088227233L, 10.4127143737989297783691126735759442322L, 0.2748768734662167365566969712751014132572L, 1.249384593879594267985829474900117552262L, 2.339628725249662110733173176958464799298L }, + { 0.13, 1.94721813462737892879360520664898861468L, 9.932951514939585997564689348668914277938L, 2.251782857706110263308481491192789580094L, 11.99260069837783378364765088349075208532L, 0.2574963206206130134352559080015524785188L, 1.250769653649645211837358227649886079707L, 2.345400881898677444465943268643293806698L }, + { 0.2, 2.142297137253033769930034117143080376893L, 13.37614224919152628643986765706565611003L, 3.535071604621394887878195503476186178003L, 45.09151212581577318204971476397941468765L, 0.1419250400262719916570182230117110839103L, 1.260560851390051221516189461818868534605L, 2.385805978441784748553427843037065714187L }, + { 0.24, 2.269461454069749995037297391358303819376L, 16.39491346167414058557668895449269328167L, 5.024462495135007919842276150900259265379L, 309.6084113035859890772945909440607406133L, 0.08069406600135331847379173643045407240867L, 1.266231911965916800406745667483278841138L, 2.408894605037846051442602142376350599479L }, +}; + +template +bool representable(long double v) +{ + return std::abs(v) <= (std::numeric_limits::max)() && (v == 0 || std::abs(v) >= (std::numeric_limits::min)()); +} + +// The condition numbers |x f'(x)/f(x)| of the row's functions, from the reference values: with +// s = 1 + xi (x - mu)/sigma and t = -logcdf, d(log t)/dx = -1/(sigma s). +struct cdf_conditions +{ + long double pdf, logpdf, cdf, ccdf, logcdf, logccdf; +}; + +inline cdf_conditions conditions(long double const (&r)[10]) +{ + long double mu = r[0], sigma = r[1], xi = r[2], x = r[3]; + long double s = 1 + xi * (x - mu) / sigma, t = -r[8], P = r[6], Q = r[7]; + long double g = std::abs(x / (sigma * s)); + cdf_conditions c; + c.pdf = g * std::abs(1 + xi - t); + c.logpdf = c.pdf / std::abs(r[5]); + c.cdf = g * t; + c.ccdf = c.cdf * P / Q; + c.logcdf = g; + c.logccdf = c.ccdf / std::abs(r[9]); + return c; +} + +template +void test_cdf() +{ + // 16 ulps times each function's condition number, which is large deep in the tails (38 for the cdf at + // xi = 0.3, x = -3, where t = 12); the condition numbers come from the long double references. + auto ulps = [](long double kappa) { return static_cast(16 * (std::max)(1.0L, kappa)); }; + for (auto const& r : cdf_data) + { + generalized_extreme_value_distribution d(static_cast(r[0]), static_cast(r[1]), static_cast(r[2])); + T x = static_cast(r[3]); + if (static_cast(x) != r[3] || static_cast(d.shape()) != r[2]) + { + continue; + } + cdf_conditions k = conditions(r); + if (representable(r[4])) { CHECK_ULP_CLOSE(static_cast(r[4]), pdf(d, x), ulps(k.pdf)); } + CHECK_ULP_CLOSE(static_cast(r[5]), logpdf(d, x), ulps(k.logpdf)); + if (representable(r[6])) { CHECK_ULP_CLOSE(static_cast(r[6]), cdf(d, x), ulps(k.cdf)); } + if (representable(r[7])) { CHECK_ULP_CLOSE(static_cast(r[7]), cdf(complement(d, x)), ulps(k.ccdf)); } + if (representable(r[8])) { CHECK_ULP_CLOSE(static_cast(r[8]), logcdf(d, x), ulps(k.logcdf)); } + if (representable(r[9])) { CHECK_ULP_CLOSE(static_cast(r[9]), logcdf(complement(d, x)), ulps(k.logccdf)); } + } +} + +template +void test_quantile() +{ + for (auto const& r : quantile_data) + { + generalized_extreme_value_distribution d(static_cast(r[0]), static_cast(r[1]), static_cast(r[2])); + T p = static_cast(r[3]); + if (static_cast(p) != r[3] || static_cast(d.shape()) != r[2]) + { + continue; + } + // mu + sigma q cancels where the quantile is near 0, multiplying the error by (|mu| + |x - mu|)/|x|: + for (int c = 0; c < 2; ++c) + { + long double x = r[4 + c]; + if (!representable(x)) + { + continue; + } + double cancellation = static_cast((std::abs(r[0]) + std::abs(x - r[0])) / std::abs(x)); + double scale = (std::max)(1.0, cancellation); +#if defined(__MINGW32__) || defined(__CYGWIN__) + // MinGW's long double pow loses about |log x| ulps for huge results: + if (std::numeric_limits::digits > 53) + { + scale = (std::max)(scale, static_cast(std::abs(std::log(std::abs(x))))); + } +#endif + unsigned ulps = static_cast(16 * scale); + CHECK_ULP_CLOSE(static_cast(x), c ? quantile(complement(d, p)) : quantile(d, p), ulps); + } + } +} + +template +void test_moments() +{ + const T eps = std::numeric_limits::epsilon(); + for (auto const& r : moment_data) + { + T xi = static_cast(r[0]); + if (static_cast(xi) != r[0]) + { + continue; + } + generalized_extreme_value_distribution d(T(0.5), T(2), xi); + CHECK_ULP_CLOSE(static_cast(r[1]), mean(d), 8); + // Gamma(1 - xi)^2 doubles tgamma's error: + CHECK_ULP_CLOSE(static_cast(r[2]), variance(d), 32); + CHECK_ULP_CLOSE(static_cast(r[2]), standard_deviation(d) * standard_deviation(d), 32); + CHECK_MOLLIFIED_CLOSE(static_cast(r[3]), skewness(d), 64 * eps); + CHECK_MOLLIFIED_CLOSE(static_cast(r[4]), kurtosis_excess(d), 512 * eps); + CHECK_MOLLIFIED_CLOSE(static_cast(r[4] + 3), kurtosis(d), 512 * eps); + CHECK_ULP_CLOSE(static_cast(r[5]), mode(d), 8); + CHECK_ULP_CLOSE(static_cast(r[6]), median(d), 8); + CHECK_ULP_CLOSE(static_cast(r[7]), entropy(d), 8); + } +} + +// The support, its endpoints, the moments' domains, and xi = 0 against extreme_value_distribution. +template +void test_edges() +{ + using boost::math::extreme_value_distribution; + const T inf = std::numeric_limits::infinity(); + generalized_extreme_value_distribution frechet(T(0.5), T(2), T(0.25)); + generalized_extreme_value_distribution weibull(T(0.5), T(2), T(-0.25)); + // mu - sigma/xi: + CHECK_EQUAL(support(frechet).first, T(-7.5)); + CHECK_EQUAL(support(weibull).second, T(8.5)); + CHECK_EQUAL(cdf(frechet, T(-7.5)), T(0)); + CHECK_EQUAL(cdf(frechet, T(-100)), T(0)); + CHECK_EQUAL(cdf(complement(frechet, T(-100))), T(1)); + CHECK_EQUAL(pdf(frechet, T(-7.5)), T(0)); + CHECK_EQUAL(logcdf(frechet, T(-100)), -inf); + CHECK_EQUAL(cdf(weibull, T(8.5)), T(1)); + CHECK_EQUAL(cdf(weibull, T(100)), T(1)); + CHECK_EQUAL(cdf(complement(weibull, T(100))), T(0)); + CHECK_EQUAL(pdf(weibull, T(8.5)), T(0)); + CHECK_EQUAL(logpdf(weibull, T(100)), -inf); + CHECK_EQUAL(quantile(frechet, T(0)), T(-7.5)); + CHECK_EQUAL(quantile(complement(weibull, T(0))), T(8.5)); + CHECK_EQUAL(cdf(frechet, inf), T(1)); + CHECK_EQUAL(cdf(weibull, -inf), T(0)); + CHECK_EQUAL(pdf(frechet, inf), T(0)); + // Just past an end point, where (x - mu)/sigma rounds onto it: + generalized_extreme_value_distribution unit_end(T(-1), T(1), T(-1)); + CHECK_EQUAL(pdf(unit_end, T(1e-17)), T(0)); + CHECK_EQUAL(logpdf(unit_end, T(1e-17)), -inf); + CHECK_EQUAL(cdf(unit_end, T(1e-17)), T(1)); + // log(1 - cdf) = -L once t = exp(-L) underflows: + generalized_extreme_value_distribution gumbel(T(0), T(1), T(0)); + if (std::numeric_limits::max_exponent10 < 300) + { + CHECK_ULP_CLOSE(T(-100), logcdf(complement(gumbel, T(100))), 1); + } + else + { + CHECK_ULP_CLOSE(T(-750), logcdf(complement(gumbel, T(750))), 1); + } + // At the upper endpoint the density is t^(1 + xi)/sigma with t = 0: + CHECK_EQUAL(pdf(generalized_extreme_value_distribution(T(0), T(2), T(-1)), T(2)), T(0.5)); + CHECK_THROW(pdf(generalized_extreme_value_distribution(T(0), T(2), T(-2)), T(1)), std::overflow_error); + // Unbounded ends: + CHECK_THROW(quantile(frechet, T(1)), std::overflow_error); + CHECK_THROW(quantile(weibull, T(0)), std::overflow_error); + // Moments that do not exist: + generalized_extreme_value_distribution heavy(T(0), T(1), T(1)); + // The mean and variance are infinite there, the skewness and kurtosis undefined: + CHECK_EQUAL(mean(heavy), inf); + CHECK_EQUAL(variance(generalized_extreme_value_distribution(T(0), T(1), T(0.5))), inf); + CHECK_EQUAL(standard_deviation(heavy), inf); + CHECK_THROW(skewness(generalized_extreme_value_distribution(T(0), T(1), T(1) / 3)), std::domain_error); + CHECK_THROW(kurtosis(generalized_extreme_value_distribution(T(0), T(1), T(0.25))), std::domain_error); + CHECK_THROW(generalized_extreme_value_distribution(T(0), T(0), T(0)), std::domain_error); + CHECK_THROW(generalized_extreme_value_distribution(T(0), T(1), inf), std::domain_error); + check_out_of_range>(T(0.5), T(2), T(0.25)); + check_out_of_range>(T(0.5), T(2), T(-0.25)); + + // xi = 0 is the Gumbel distribution, and tiny xi tends to it; two implementations, deep in the tail at x = 30: + generalized_extreme_value_distribution g(T(0.5), T(2), T(0)); + extreme_value_distribution e(T(0.5), T(2)); + for (T x : { T(-3), T(0.5), T(4), T(30) }) + { + CHECK_ULP_CLOSE(pdf(e, x), pdf(g, x), 16); + CHECK_ULP_CLOSE(cdf(e, x), cdf(g, x), 4); + CHECK_ULP_CLOSE(cdf(complement(e, x)), cdf(complement(g, x)), 4); + } + for (T p : { T(0.001), T(0.5), T(0.999) }) + { + CHECK_ULP_CLOSE(quantile(e, p), quantile(g, p), 4); + } + CHECK_ULP_CLOSE(skewness(e), skewness(g), 0); + CHECK_ULP_CLOSE(kurtosis(e), kurtosis(g), 0); + generalized_extreme_value_distribution tiny(T(0.5), T(2), std::numeric_limits::denorm_min()); + CHECK_ULP_CLOSE(mean(e), mean(tiny), 1); + CHECK_ULP_CLOSE(cdf(e, T(4)), cdf(tiny, T(4)), 4); + CHECK_ULP_CLOSE(quantile(e, T(0.25)), quantile(tiny, T(0.25)), 4); +} + +#ifdef BOOST_MATH_TEST_GEV_MULTIPRECISION +// The moment series in many digits, against the closed forms in more, and the quantile's round trip +// (away from p -> 0, where cdf(quantile(p)) multiplies the error by -log p). +void test_multiprecision() +{ + using T = boost::multiprecision::cpp_bin_float_50; + using U = boost::multiprecision::cpp_bin_float_100; + for (const char* s : { "-0.1", "-1e-10", "1e-10", "0.0625", "0.15", "0.2" }) + { + T xi(s); + U k(xi); + generalized_extreme_value_distribution d(T(0), T(1), xi); + U g1 = boost::math::tgamma(1 - k), g2 = boost::math::tgamma(1 - 2 * k), g3 = boost::math::tgamma(1 - 3 * k), g4 = boost::math::tgamma(1 - 4 * k); + U var = (g2 - g1 * g1) / (k * k); + U skew = (k > 0 ? 1 : -1) * (g3 - 3 * g1 * g2 + 2 * g1 * g1 * g1) / pow(g2 - g1 * g1, U(1.5)); + U kex = (g4 - 4 * g1 * g3 + 6 * g2 * g1 * g1 - 3 * g1 * g1 * g1 * g1) / ((g2 - g1 * g1) * (g2 - g1 * g1)) - 3; + // tgamma itself is ~20 ulps at 50 digits, and the variance has Gamma(1 - xi)^2: + CHECK_ULP_CLOSE(T(var), variance(d), 128); + CHECK_MOLLIFIED_CLOSE(T(skew), skewness(d), T(64 * std::numeric_limits::epsilon())); + CHECK_MOLLIFIED_CLOSE(T(kex), kurtosis_excess(d), T(512 * std::numeric_limits::epsilon())); + for (T p : { T(0.001), T(0.3), T(0.9) }) + { + CHECK_ULP_CLOSE(p, cdf(d, quantile(d, p)), 64); + } + } +} +#endif + +int main() +{ + test_cdf(); + test_cdf(); + test_cdf(); + test_quantile(); + test_quantile(); + test_quantile(); + test_moments(); + test_moments(); + test_moments(); + test_edges(); + test_edges(); + test_edges(); +#ifdef BOOST_MATH_TEST_GEV_MULTIPRECISION + test_multiprecision(); +#endif + return boost::math::test::report_errors(); +} diff --git a/test/test_negative_binomial.cpp b/test/test_negative_binomial.cpp index 271e4f0cee..827dff4da5 100644 --- a/test/test_negative_binomial.cpp +++ b/test/test_negative_binomial.cpp @@ -31,32 +31,25 @@ #endif #include -#include "../include_private/boost/math/tools/test.hpp" +#include "math_unit_test.hpp" #ifndef BOOST_MATH_NO_REAL_CONCEPT_TESTS #include // for real_concept -using ::boost::math::concepts::real_concept; #endif -#include // for negative_binomial_distribution -using boost::math::negative_binomial_distribution; - -#include - using boost::math::lgamma; // log gamma - -#define BOOST_TEST_MAIN -#include // for test_main -#include // for BOOST_CHECK_CLOSE +#include +#include #include "table_type.hpp" #include "test_out_of_range.hpp" -#include -using std::cout; -using std::endl; -using std::setprecision; -using std::showpoint; +#include #include -using std::numeric_limits; +#include + +using boost::math::negative_binomial_distribution; + +// The Boost.Test checks these replaced took tolerances in percent; the tolerances below are +// fractions, the same values divided by 100. template void test_spot( // Test a single spot value against 'known good' values. @@ -67,124 +60,50 @@ void test_spot( // Test a single spot value against 'known good' values. RealType Q, // Complement of CDF. RealType tol) // Test tolerance. { - boost::math::negative_binomial_distribution bn(N, p); - BOOST_CHECK_EQUAL(N, bn.successes()); - BOOST_CHECK_EQUAL(p, bn.success_fraction()); - BOOST_CHECK_CLOSE( - cdf(bn, k), P, tol); + negative_binomial_distribution bn(N, p); + CHECK_EQUAL(N, bn.successes()); + CHECK_EQUAL(p, bn.success_fraction()); + CHECK_RELATIVE_ERROR(P, cdf(bn, k), tol); if((P < 0.99) && (Q < 0.99)) { // We can only check this if P is not too close to 1, // so that we can guarantee that Q is free of error: - // - BOOST_CHECK_CLOSE( - cdf(complement(bn, k)), Q, tol); - if(k != 0) - { - BOOST_CHECK_CLOSE( - quantile(bn, P), k, tol); - } - else - { - // Just check quantile is very small: - if((std::numeric_limits::max_exponent <= std::numeric_limits::max_exponent) - && (boost::is_floating_point::value)) - { - // Limit where this is checked: if exponent range is very large we may - // run out of iterations in our root finding algorithm. - BOOST_CHECK(quantile(bn, P) < boost::math::tools::epsilon() * 10); - } - } + CHECK_RELATIVE_ERROR(Q, cdf(complement(bn, k)), tol); + // Limit where a zero quantile is checked: if exponent range is very large we may + // run out of iterations in our root finding algorithm. + bool check_zero = (std::numeric_limits::max_exponent <= std::numeric_limits::max_exponent) + && std::is_floating_point::value; if(k != 0) { - BOOST_CHECK_CLOSE( - quantile(complement(bn, Q)), k, tol); + CHECK_RELATIVE_ERROR(k, quantile(bn, P), tol); + CHECK_RELATIVE_ERROR(k, quantile(complement(bn, Q)), tol); } - else + else if(check_zero) { // Just check quantile is very small: - if((std::numeric_limits::max_exponent <= std::numeric_limits::max_exponent) - && (boost::is_floating_point::value)) - { - // Limit where this is checked: if exponent range is very large we may - // run out of iterations in our root finding algorithm. - BOOST_CHECK(quantile(complement(bn, Q)) < boost::math::tools::epsilon() * 10); - } + CHECK_LE(quantile(bn, P), boost::math::tools::epsilon() * 10); + CHECK_LE(quantile(complement(bn, Q)), boost::math::tools::epsilon() * 10); } // estimate success ratio: - BOOST_CHECK_CLOSE( - negative_binomial_distribution::find_lower_bound_on_p( - N+k, N, P), - p, tol); + CHECK_RELATIVE_ERROR(p, negative_binomial_distribution::find_lower_bound_on_p(N+k, N, P), tol); // Note we bump up the sample size here, purely for the sake of the test, // internally the function has to adjust the sample size so that we get // the right upper bound, our test undoes this, so we can verify the result. - BOOST_CHECK_CLOSE( - negative_binomial_distribution::find_upper_bound_on_p( - N+k+1, N, Q), - p, tol); + CHECK_RELATIVE_ERROR(p, negative_binomial_distribution::find_upper_bound_on_p(N+k+1, N, Q), tol); - if(Q < P) - { - // - // We check two things here, that the upper and lower bounds - // are the right way around, and that they do actually bracket - // the naive estimate of p = successes / (sample size) - // - BOOST_CHECK( - negative_binomial_distribution::find_lower_bound_on_p( - N+k, N, Q) - <= - negative_binomial_distribution::find_upper_bound_on_p( - N+k, N, Q) - ); - BOOST_CHECK( - negative_binomial_distribution::find_lower_bound_on_p( - N+k, N, Q) - <= - N / (N+k) - ); - BOOST_CHECK( - N / (N+k) - <= - negative_binomial_distribution::find_upper_bound_on_p( - N+k, N, Q) - ); - } - else - { - // As above but when P is small. - BOOST_CHECK( - negative_binomial_distribution::find_lower_bound_on_p( - N+k, N, P) - <= - negative_binomial_distribution::find_upper_bound_on_p( - N+k, N, P) - ); - BOOST_CHECK( - negative_binomial_distribution::find_lower_bound_on_p( - N+k, N, P) - <= - N / (N+k) - ); - BOOST_CHECK( - N / (N+k) - <= - negative_binomial_distribution::find_upper_bound_on_p( - N+k, N, P) - ); - } + // Check that the upper and lower bounds are the right way around, and that they + // bracket the naive estimate of p = successes / (sample size); with Q when Q < P, else with P: + RealType alpha = (Q < P) ? Q : P; + RealType lower = negative_binomial_distribution::find_lower_bound_on_p(N+k, N, alpha); + RealType upper = negative_binomial_distribution::find_upper_bound_on_p(N+k, N, alpha); + CHECK_LE(lower, upper); + CHECK_LE(lower, RealType(N / (N+k))); + CHECK_LE(RealType(N / (N+k)), upper); // Estimate sample size: - BOOST_CHECK_CLOSE( - negative_binomial_distribution::find_minimum_number_of_trials( - k, p, P), - N+k, tol); - BOOST_CHECK_CLOSE( - negative_binomial_distribution::find_maximum_number_of_trials( - k, p, Q), - N+k, tol); + CHECK_RELATIVE_ERROR(N+k, negative_binomial_distribution::find_minimum_number_of_trials(k, p, P), tol); + CHECK_RELATIVE_ERROR(N+k, negative_binomial_distribution::find_maximum_number_of_trials(k, p, Q), tol); // Double check consistency of CDF and PDF by computing the finite sum: RealType sum = 0; @@ -192,31 +111,45 @@ void test_spot( // Test a single spot value against 'known good' values. { sum += pdf(bn, RealType(i)); } - BOOST_CHECK_CLOSE(sum, P, tol); - + CHECK_RELATIVE_ERROR(P, sum, tol); // Complement is not possible since sum is to infinity. - } // + } } // test_spot +// The relative condition number |P / (x dP/dx)| of the quantile x of P, from a central difference +// of the cdf: the quantile's relative error is the cdf's times this, which is large near x = 0. +template +RealType quantile_condition_number(Dist const& d, RealType x, RealType P, bool complemented) +{ + using std::abs; + using std::pow; + if(x == 0) + { + return 1; + } + // Not tools::cbrt_epsilon, whose function-local static a SYCL kernel can't use: + RealType h = x * pow(boost::math::tools::epsilon(), RealType(1) / 3); + RealType up = complemented ? cdf(complement(d, RealType(x + h))) : cdf(d, RealType(x + h)); + RealType down = complemented ? cdf(complement(d, RealType(x - h))) : cdf(d, RealType(x - h)); + RealType kappa = abs(P * 2 * h / (x * (up - down))); + return kappa > 1 ? kappa : RealType(1); +} + template // Any floating-point type RealType. void test_spots(RealType) { // Basic sanity checks, test data is to double precision only - // so set tolerance to 1000 eps expressed as a percent, or - // 1000 eps of type double expressed as a percent, whichever - // is the larger. - + // so set tolerance to 100000 eps, or 100000 eps of type double, + // whichever is the larger. RealType tolerance = (std::max) (boost::math::tools::epsilon(), static_cast(std::numeric_limits::epsilon())); - tolerance *= 100 * 100000.0f; + tolerance *= 100000.0f; + // BOOST_CHECK_SMALL took the percentage tolerance above as an absolute bound: + RealType small_bound = tolerance * 100; - cout << "Tolerance = " << tolerance << "%." << endl; - - RealType tol1eps = boost::math::tools::epsilon() * 2; // Very tight, suit exact values. - //RealType tol2eps = boost::math::tools::epsilon() * 2; // Tight, suit exact values. - RealType tol5eps = boost::math::tools::epsilon() * 5; // Wider 5 epsilon. - cout << "Tolerance 5 eps = " << tol5eps << "%." << endl; + RealType tol1eps = boost::math::tools::epsilon() * 2 / 100; // Very tight, suit exact values. + RealType tol5eps = boost::math::tools::epsilon() * 5 / 100; // Wider 5 epsilon. // Sources of spot test values: @@ -229,521 +162,174 @@ void test_spots(RealType) // And functions.wolfram.com - using boost::math::negative_binomial_distribution; - using ::boost::math::negative_binomial; - using ::boost::math::cdf; - using ::boost::math::pdf; + using ::boost::math::cdf; + using ::boost::math::pdf; // Test negative binomial using cdf spot values from MathCAD cdf = pnbinom(k, r, p). // These test quantiles and complements as well. + // Arguments: successes r, failures k, success fraction p, CDF P, complement Q = 1 - P. test_spot( // pnbinom(1,2,0.5) = 0.5 - static_cast(2), // successes r - static_cast(1), // Number of failures, k - static_cast(0.5), // Probability of success as fraction, p - static_cast(0.5), // Probability of result (CDF), P - static_cast(0.5), // complement CCDF Q = 1 - P - tolerance); + static_cast(2), static_cast(1), static_cast(0.5), + static_cast(0.5), static_cast(0.5), tolerance); test_spot( // pbinom(0, 2, 0.25) - static_cast(2), // successes r - static_cast(0), // Number of failures, k - static_cast(0.25), - static_cast(0.0625), // Probability of result (CDF), P - static_cast(0.9375), // Q = 1 - P - tolerance); + static_cast(2), static_cast(0), static_cast(0.25), + static_cast(0.0625), static_cast(0.9375), tolerance); test_spot( // pbinom(48,8,0.25) - static_cast(8), // successes r - static_cast(48), // Number of failures, k - static_cast(0.25), // Probability of success, p - static_cast(9.826582228110670E-1), // Probability of result (CDF), P - static_cast(1 - 9.826582228110670E-1), // Q = 1 - P - tolerance); + static_cast(8), static_cast(48), static_cast(0.25), + static_cast(9.826582228110670E-1), static_cast(1 - 9.826582228110670E-1), tolerance); test_spot( // pbinom(2,5,0.4) - static_cast(5), // successes r - static_cast(2), // Number of failures, k - static_cast(0.4), // Probability of success, p - static_cast(9.625600000000020E-2), // Probability of result (CDF), P - static_cast(1 - 9.625600000000020E-2), // Q = 1 - P - tolerance); + static_cast(5), static_cast(2), static_cast(0.4), + static_cast(9.625600000000020E-2), static_cast(1 - 9.625600000000020E-2), tolerance); test_spot( // pbinom(10,100,0.9) - static_cast(100), // successes r - static_cast(10), // Number of failures, k - static_cast(0.9), // Probability of success, p - static_cast(4.535522887695670E-1), // Probability of result (CDF), P - static_cast(1 - 4.535522887695670E-1), // Q = 1 - P - tolerance); + static_cast(100), static_cast(10), static_cast(0.9), + static_cast(4.535522887695670E-1), static_cast(1 - 4.535522887695670E-1), tolerance); test_spot( // pbinom(1,100,0.991) - static_cast(100), // successes r - static_cast(1), // Number of failures, k - static_cast(0.991), // Probability of success, p - static_cast(7.693413044217000E-1), // Probability of result (CDF), P - static_cast(1 - 7.693413044217000E-1), // Q = 1 - P - tolerance); + static_cast(100), static_cast(1), static_cast(0.991), + static_cast(7.693413044217000E-1), static_cast(1 - 7.693413044217000E-1), tolerance); test_spot( // pbinom(10,100,0.991) - static_cast(100), // successes r - static_cast(10), // Number of failures, k - static_cast(0.991), // Probability of success, p - static_cast(9.999999940939000E-1), // Probability of result (CDF), P - static_cast(1 - 9.999999940939000E-1), // Q = 1 - P - tolerance); - -if(std::numeric_limits::is_specialized) -{ // An extreme value test that takes 3 minutes using the real concept type - // for which numeric_limits::is_specialized == false, deliberately - // and for which there is no Lanczos approximation defined (also deliberately) - // giving a very slow computation, but with acceptable accuracy. - // A possible enhancement might be to use a normal approximation for - // extreme values, but this is not implemented. - test_spot( // pbinom(100000,100,0.001) - static_cast(100), // successes r - static_cast(100000), // Number of failures, k - static_cast(0.001), // Probability of success, p - static_cast(5.173047534260320E-1), // Probability of result (CDF), P - static_cast(1 - 5.173047534260320E-1), // Q = 1 - P - tolerance*1000); // *1000 is OK 0.51730475350664229 versus - - // functions.wolfram.com - // for I[0.001](100, 100000+1) gives: - // Wolfram 0.517304753506834882009032744488738352004003696396461766326713 - // JM nonLanczos 0.51730475350664229 differs at the 13th decimal digit. - // MathCAD 0.51730475342603199 differs at 10th decimal digit. - - // Error tests: - check_out_of_range >(20, 0.5); - BOOST_MATH_CHECK_THROW(negative_binomial_distribution(0, 0.5), std::domain_error); - BOOST_MATH_CHECK_THROW(negative_binomial_distribution(-2, 0.5), std::domain_error); - BOOST_MATH_CHECK_THROW(negative_binomial_distribution(20, -0.5), std::domain_error); - BOOST_MATH_CHECK_THROW(negative_binomial_distribution(20, 1.5), std::domain_error); -} - // End of single spot tests using RealType - - - // Tests on PDF: - BOOST_CHECK_CLOSE( - pdf(negative_binomial_distribution(static_cast(2), static_cast(0.5)), - static_cast(0) ), // k = 0. - static_cast(0.25), // 0 - tolerance); - - BOOST_CHECK_CLOSE( - pdf(negative_binomial_distribution(static_cast(4), static_cast(0.5)), - static_cast(0)), // k = 0. - static_cast(0.0625), // exact 1/16 - tolerance); - - BOOST_CHECK_CLOSE( - pdf(negative_binomial_distribution(static_cast(20), static_cast(0.25)), - static_cast(0)), // k = 0 - static_cast(9.094947017729270E-13), // pbinom(0,20,0.25) = 9.094947017729270E-13 - tolerance); - - BOOST_CHECK_CLOSE( - pdf(negative_binomial_distribution(static_cast(20), static_cast(0.2)), - static_cast(0)), // k = 0 - static_cast(1.0485760000000003e-014), // MathCAD 1.048576000000000E-14 - tolerance); - - BOOST_CHECK_CLOSE( - pdf(negative_binomial_distribution(static_cast(10), static_cast(0.1)), - static_cast(0)), // k = 0. - static_cast(1e-10), // MathCAD says zero, but suffers cancellation error? - tolerance); - - BOOST_CHECK_CLOSE( - pdf(negative_binomial_distribution(static_cast(20), static_cast(0.1)), - static_cast(0)), // k = 0. - static_cast(1e-20), // MathCAD says zero, but suffers cancellation error? - tolerance); - - - BOOST_CHECK_CLOSE( // . - pdf(negative_binomial_distribution(static_cast(20), static_cast(0.9)), - static_cast(0)), // k. - static_cast(1.215766545905690E-1), // k=20 p = 0.9 - tolerance); + static_cast(100), static_cast(10), static_cast(0.991), + static_cast(9.999999940939000E-1), static_cast(1 - 9.999999940939000E-1), tolerance); + + if(std::numeric_limits::is_specialized) + { // An extreme value test that takes 3 minutes using the real concept type + // for which numeric_limits::is_specialized == false, deliberately + // and for which there is no Lanczos approximation defined (also deliberately) + // giving a very slow computation, but with acceptable accuracy. + // A possible enhancement might be to use a normal approximation for + // extreme values, but this is not implemented. + test_spot( // pbinom(100000,100,0.001) + static_cast(100), static_cast(100000), static_cast(0.001), + static_cast(5.173047534260320E-1), static_cast(1 - 5.173047534260320E-1), + tolerance*1000); // *1000 is OK 0.51730475350664229 versus + + // functions.wolfram.com + // for I[0.001](100, 100000+1) gives: + // Wolfram 0.517304753506834882009032744488738352004003696396461766326713 + // JM nonLanczos 0.51730475350664229 differs at the 13th decimal digit. + // MathCAD 0.51730475342603199 differs at 10th decimal digit. + + // Error tests: + check_out_of_range >(20, 0.5); + CHECK_THROW(negative_binomial_distribution(0, 0.5), std::domain_error); + CHECK_THROW(negative_binomial_distribution(-2, 0.5), std::domain_error); + CHECK_THROW(negative_binomial_distribution(20, -0.5), std::domain_error); + CHECK_THROW(negative_binomial_distribution(20, 1.5), std::domain_error); + } + // End of single spot tests using RealType + + auto nb = [](double r, double p) { return negative_binomial_distribution(static_cast(r), static_cast(p)); }; + + // Tests on PDF at k = 0: + CHECK_RELATIVE_ERROR(static_cast(0.25), pdf(nb(2, 0.5), static_cast(0)), tolerance); + CHECK_RELATIVE_ERROR(static_cast(0.0625), pdf(nb(4, 0.5), static_cast(0)), tolerance); // exact 1/16 + CHECK_RELATIVE_ERROR(static_cast(9.094947017729270E-13), pdf(nb(20, 0.25), static_cast(0)), tolerance); // pbinom(0,20,0.25) + CHECK_RELATIVE_ERROR(static_cast(1.0485760000000003e-014), pdf(nb(20, 0.2), static_cast(0)), tolerance); // MathCAD 1.048576000000000E-14 + CHECK_RELATIVE_ERROR(static_cast(1e-10), pdf(nb(10, 0.1), static_cast(0)), tolerance); // MathCAD says zero, but suffers cancellation error? + CHECK_RELATIVE_ERROR(static_cast(1e-20), pdf(nb(20, 0.1), static_cast(0)), tolerance); // MathCAD says zero, but suffers cancellation error? + CHECK_RELATIVE_ERROR(static_cast(1.215766545905690E-1), pdf(nb(20, 0.9), static_cast(0)), tolerance); // Tests on cdf: // MathCAD pbinom k, r, p) == failures, successes, probability. - - BOOST_CHECK_CLOSE(cdf( - negative_binomial_distribution(static_cast(2), static_cast(0.5)), // successes = 2,prob 0.25 - static_cast(0) ), // k = 0 - static_cast(0.25), // probability 1/4 - tolerance); - - BOOST_CHECK_CLOSE(cdf(complement( - negative_binomial_distribution(static_cast(2), static_cast(0.5)), // successes = 2,prob 0.25 - static_cast(0) )), // k = 0 - static_cast(0.75), // probability 3/4 - tolerance); - BOOST_CHECK_CLOSE( // k = 1. - cdf(negative_binomial_distribution(static_cast(20), static_cast(0.25)), - static_cast(1)), // k =1. - static_cast(1.455191522836700E-11), - tolerance); - - BOOST_CHECK_SMALL( // Check within an epsilon with CHECK_SMALL - cdf(negative_binomial_distribution(static_cast(20), static_cast(0.25)), - static_cast(1)) - - static_cast(1.455191522836700E-11), - tolerance ); + CHECK_RELATIVE_ERROR(static_cast(0.25), cdf(nb(2, 0.5), static_cast(0)), tolerance); + CHECK_RELATIVE_ERROR(static_cast(0.75), cdf(complement(nb(2, 0.5), static_cast(0))), tolerance); + CHECK_RELATIVE_ERROR(static_cast(1.455191522836700E-11), cdf(nb(20, 0.25), static_cast(1)), tolerance); + CHECK_ABSOLUTE_ERROR(static_cast(1.455191522836700E-11), cdf(nb(20, 0.25), static_cast(1)), small_bound); // Some exact (probably - judging by trailing zeros) values. - BOOST_CHECK_CLOSE( - cdf(negative_binomial_distribution(static_cast(8), static_cast(0.25)), - static_cast(0)), // k. - static_cast(1.525878906250000E-5), - tolerance); - - BOOST_CHECK_CLOSE( - cdf(negative_binomial_distribution(static_cast(8), static_cast(0.25)), - static_cast(0)), // k. - static_cast(1.525878906250000E-5), - tolerance); - - BOOST_CHECK_SMALL( - cdf(negative_binomial_distribution(static_cast(8), static_cast(0.25)), - static_cast(0)) - - static_cast(1.525878906250000E-5), - tolerance ); - - BOOST_CHECK_CLOSE( // k = 1. - cdf(negative_binomial_distribution(static_cast(8), static_cast(0.25)), - static_cast(1)), // k. - static_cast(1.068115234375010E-4), - tolerance); - - BOOST_CHECK_CLOSE( // k = 2. - cdf(negative_binomial_distribution(static_cast(8), static_cast(0.25)), - static_cast(2)), // k. - static_cast(4.158020019531300E-4), - tolerance); - - BOOST_CHECK_CLOSE( // k = 3. - cdf(negative_binomial_distribution(static_cast(8), static_cast(0.25)), - static_cast(3)), // k.bristow - static_cast(1.188278198242200E-3), - tolerance); - - BOOST_CHECK_CLOSE( // k = 4. - cdf(negative_binomial_distribution(static_cast(8), static_cast(0.25)), - static_cast(4)), // k. - static_cast(2.781510353088410E-3), - tolerance); - - BOOST_CHECK_CLOSE( // k = 5. - cdf(negative_binomial_distribution(static_cast(8), static_cast(0.25)), - static_cast(5)), // k. - static_cast(5.649328231811500E-3), - tolerance); - - BOOST_CHECK_CLOSE( // k = 6. - cdf(negative_binomial_distribution(static_cast(8), static_cast(0.25)), - static_cast(6)), // k. - static_cast(1.030953228473680E-2), - tolerance); - - BOOST_CHECK_CLOSE( // k = 7. - cdf(negative_binomial_distribution(static_cast(8), static_cast(0.25)), - static_cast(7)), // k. - static_cast(1.729983836412430E-2), - tolerance); - - BOOST_CHECK_CLOSE( // k = 8. - cdf(negative_binomial_distribution(static_cast(8), static_cast(0.25)), - static_cast(8)), // k = n. - static_cast(2.712995628826370E-2), - tolerance); - - BOOST_CHECK_CLOSE( // - cdf(negative_binomial_distribution(static_cast(8), static_cast(0.25)), - static_cast(48)), // k - static_cast(9.826582228110670E-1), - tolerance); - - BOOST_CHECK_CLOSE( // - cdf(negative_binomial_distribution(static_cast(8), static_cast(0.25)), - static_cast(64)), // k - static_cast(9.990295004935590E-1), - tolerance); - - BOOST_CHECK_CLOSE( // - cdf(negative_binomial_distribution(static_cast(5), static_cast(0.4)), - static_cast(26)), // k - static_cast(9.989686246611190E-1), - tolerance); - - BOOST_CHECK_CLOSE( // - cdf(negative_binomial_distribution(static_cast(5), static_cast(0.4)), - static_cast(2)), // k failures - static_cast(9.625600000000020E-2), - tolerance); - - BOOST_CHECK_CLOSE( // - cdf(negative_binomial_distribution(static_cast(50), static_cast(0.9)), - static_cast(20)), // k - static_cast(9.999970854144170E-1), - tolerance); - - BOOST_CHECK_CLOSE( // - cdf(negative_binomial_distribution(static_cast(500), static_cast(0.7)), - static_cast(200)), // k - static_cast(2.172846379930550E-1), - tolerance* 2); - - BOOST_CHECK_CLOSE( // - cdf(negative_binomial_distribution(static_cast(50), static_cast(0.7)), - static_cast(20)), // k - static_cast(4.550203671301790E-1), - tolerance); + CHECK_RELATIVE_ERROR(static_cast(1.525878906250000E-5), cdf(nb(8, 0.25), static_cast(0)), tolerance); + CHECK_ABSOLUTE_ERROR(static_cast(1.525878906250000E-5), cdf(nb(8, 0.25), static_cast(0)), small_bound); + const double eight_quarter[][2] = { + {1, 1.068115234375010E-4}, {2, 4.158020019531300E-4}, {3, 1.188278198242200E-3}, {4, 2.781510353088410E-3}, + {5, 5.649328231811500E-3}, {6, 1.030953228473680E-2}, {7, 1.729983836412430E-2}, {8, 2.712995628826370E-2}, + {48, 9.826582228110670E-1}, {64, 9.990295004935590E-1} }; + for(auto const& row : eight_quarter) + { + CHECK_RELATIVE_ERROR(static_cast(row[1]), cdf(nb(8, 0.25), static_cast(row[0])), tolerance); + } + CHECK_RELATIVE_ERROR(static_cast(9.989686246611190E-1), cdf(nb(5, 0.4), static_cast(26)), tolerance); + CHECK_RELATIVE_ERROR(static_cast(9.625600000000020E-2), cdf(nb(5, 0.4), static_cast(2)), tolerance); + CHECK_RELATIVE_ERROR(static_cast(9.999970854144170E-1), cdf(nb(50, 0.9), static_cast(20)), tolerance); + CHECK_RELATIVE_ERROR(static_cast(2.172846379930550E-1), cdf(nb(500, 0.7), static_cast(200)), tolerance * 2); + CHECK_RELATIVE_ERROR(static_cast(4.550203671301790E-1), cdf(nb(50, 0.7), static_cast(20)), tolerance); // Tests of other functions, mean and other moments ... - negative_binomial_distribution dist(static_cast(8), static_cast(0.25)); + negative_binomial_distribution dist = nb(8, 0.25); using namespace std; // ADL of std names. - // mean: - BOOST_CHECK_CLOSE( - mean(dist), static_cast(8 * (1 - 0.25) /0.25), tol5eps); - BOOST_CHECK_CLOSE( - mode(dist), static_cast(21), tol1eps); - // variance: - BOOST_CHECK_CLOSE( - variance(dist), static_cast(8 * (1 - 0.25) / (0.25 * 0.25)), tol5eps); - // std deviation: - BOOST_CHECK_CLOSE( - standard_deviation(dist), // 9.79795897113271239270 - static_cast(9.797958971132712392789136298823565567864L), // using functions.wolfram.com - // 9.79795897113271152534 == sqrt(8 * (1 - 0.25) / (0.25 * 0.25))) - tol5eps * 100); - BOOST_CHECK_CLOSE( - skewness(dist), // - static_cast(0.71443450831176036), - // using http://mathworld.wolfram.com/skewness.html - tolerance); - BOOST_CHECK_CLOSE( - kurtosis_excess(dist), // - static_cast(0.7604166666666666666666666666666666666666L), // using Wikipedia Kurtosis(excess) formula - tol5eps * 100); - BOOST_CHECK_CLOSE( - kurtosis(dist), // true - static_cast(3.76041666666666666666666666666666666666666L), // - tol5eps * 100); - // hazard: + CHECK_RELATIVE_ERROR(static_cast(8 * (1 - 0.25) /0.25), mean(dist), tol5eps); + CHECK_RELATIVE_ERROR(static_cast(21), mode(dist), tol1eps); + CHECK_RELATIVE_ERROR(static_cast(8 * (1 - 0.25) / (0.25 * 0.25)), variance(dist), tol5eps); + CHECK_RELATIVE_ERROR(static_cast(9.797958971132712392789136298823565567864L), standard_deviation(dist), tol5eps * 100); // using functions.wolfram.com + CHECK_RELATIVE_ERROR(static_cast(0.71443450831176036), skewness(dist), tolerance); // using http://mathworld.wolfram.com/skewness.html + CHECK_RELATIVE_ERROR(static_cast(0.7604166666666666666666666666666666666666L), kurtosis_excess(dist), tol5eps * 100); // using Wikipedia Kurtosis(excess) formula + CHECK_RELATIVE_ERROR(static_cast(3.76041666666666666666666666666666666666666L), kurtosis(dist), tol5eps * 100); RealType x = static_cast(0.125); - BOOST_CHECK_CLOSE( - hazard(dist, x) - , pdf(dist, x) / cdf(complement(dist, x)), tol5eps); - // cumulative hazard: - BOOST_CHECK_CLOSE( - chf(dist, x), -log(cdf(complement(dist, x))), tol5eps); - // coefficient_of_variation: - BOOST_CHECK_CLOSE( - coefficient_of_variation(dist) - , standard_deviation(dist) / mean(dist), tol5eps); + CHECK_RELATIVE_ERROR(pdf(dist, x) / cdf(complement(dist, x)), hazard(dist, x), tol5eps); + CHECK_RELATIVE_ERROR(-log(cdf(complement(dist, x))), chf(dist, x), tol5eps); + CHECK_RELATIVE_ERROR(standard_deviation(dist) / mean(dist), coefficient_of_variation(dist), tol5eps); // Special cases for PDF: - BOOST_CHECK_EQUAL( - pdf( - negative_binomial_distribution(static_cast(8), static_cast(0)), // - static_cast(0)), - static_cast(0) ); - - BOOST_CHECK_EQUAL( - pdf( - negative_binomial_distribution(static_cast(8), static_cast(0)), - static_cast(0.0001)), - static_cast(0) ); - - BOOST_CHECK_EQUAL( - pdf( - negative_binomial_distribution(static_cast(8), static_cast(1)), - static_cast(0.001)), - static_cast(0) ); - - BOOST_CHECK_EQUAL( - pdf( - negative_binomial_distribution(static_cast(8), static_cast(1)), - static_cast(8)), - static_cast(0) ); - - BOOST_CHECK_SMALL( - pdf( - negative_binomial_distribution(static_cast(2), static_cast(0.25)), - static_cast(0))- - static_cast(0.0625), - 2 * boost::math::tools::epsilon() ); // Expect exact, but not quite. - // numeric_limits::epsilon()); // Not suitable for real concept! + CHECK_EQUAL(pdf(nb(8, 0), static_cast(0)), static_cast(0)); + CHECK_EQUAL(pdf(nb(8, 0), static_cast(0.0001)), static_cast(0)); + CHECK_EQUAL(pdf(nb(8, 1), static_cast(0.001)), static_cast(0)); + CHECK_EQUAL(pdf(nb(8, 1), static_cast(8)), static_cast(0)); + // Expect exact, but not quite: + CHECK_ABSOLUTE_ERROR(static_cast(0.0625), pdf(nb(2, 0.25), static_cast(0)), RealType(2 * boost::math::tools::epsilon())); // Quantile boundary cases checks: - BOOST_CHECK_EQUAL( - quantile( // zero P < cdf(0) so should be exactly zero. - negative_binomial_distribution(static_cast(8), static_cast(0.25)), - static_cast(0)), - static_cast(0)); - - BOOST_CHECK_EQUAL( - quantile( // min P < cdf(0) so should be exactly zero. - negative_binomial_distribution(static_cast(8), static_cast(0.25)), - static_cast(boost::math::tools::min_value())), - static_cast(0)); - - BOOST_CHECK_CLOSE_FRACTION( - quantile( // Small P < cdf(0) so should be near zero. - negative_binomial_distribution(static_cast(8), static_cast(0.25)), - static_cast(boost::math::tools::epsilon())), // - static_cast(0), - tol5eps); - - BOOST_CHECK_CLOSE( - quantile( // Small P < cdf(0) so should be exactly zero. - negative_binomial_distribution(static_cast(8), static_cast(0.25)), - static_cast(0.0001)), - static_cast(0.95854156929288470), - tolerance); - - //BOOST_CHECK( // Fails with overflow for real_concept - //quantile( // Small P near 1 so k failures should be big. - //negative_binomial_distribution(static_cast(8), static_cast(0.25)), - //static_cast(1 - boost::math::tools::epsilon())) <= - //static_cast(189.56999032670058) // 106.462769 for float - //); + CHECK_EQUAL(quantile(dist, static_cast(0)), static_cast(0)); // zero P < cdf(0) so should be exactly zero. + CHECK_EQUAL(quantile(dist, boost::math::tools::min_value()), static_cast(0)); // min P < cdf(0) + // Small P < cdf(0) so should be near zero; a relative check against zero requires zero: + CHECK_RELATIVE_ERROR(static_cast(0), quantile(dist, boost::math::tools::epsilon()), tol5eps); + CHECK_RELATIVE_ERROR(static_cast(0.95854156929288470), quantile(dist, static_cast(0.0001)), tolerance); if(std::numeric_limits::has_infinity) - { // BOOST_CHECK tests for infinity using std::numeric_limits<>::infinity() - // Note that infinity is not implemented for real_concept, so these tests - // are only done for types, like built-in float, double.. that have infinity. - // Note that these assume that BOOST_MATH_OVERFLOW_ERROR_POLICY is NOT throw_on_error. - // #define BOOST_MATH_THROW_ON_OVERFLOW_POLICY == throw_on_error would throw here. - // #define BOOST_MAT_DOMAIN_ERROR_POLICY IS defined throw_on_error, - // so the throw path of error handling is tested below with BOOST_MATH_CHECK_THROW tests. - - BOOST_CHECK( - quantile( // At P == 1 so k failures should be infinite. - negative_binomial_distribution(static_cast(8), static_cast(0.25)), - static_cast(1)) == - //static_cast(boost::math::tools::infinity()) - static_cast(std::numeric_limits::infinity()) ); - - BOOST_CHECK_EQUAL( - quantile( // At 1 == P so should be infinite. - negative_binomial_distribution(static_cast(8), static_cast(0.25)), - static_cast(1)), // - std::numeric_limits::infinity() ); - - BOOST_CHECK_EQUAL( - quantile(complement( // Q zero 1 so P == 1 < cdf(0) so should be exactly infinity. - negative_binomial_distribution(static_cast(8), static_cast(0.25)), - static_cast(0))), - std::numeric_limits::infinity() ); - } // test for infinity using std::numeric_limits<>::infinity() + { // Infinity is not implemented for real_concept, so these tests are only done for built-in types. + // These assume that BOOST_MATH_OVERFLOW_ERROR_POLICY is NOT throw_on_error. + CHECK_EQUAL(quantile(dist, static_cast(1)), std::numeric_limits::infinity()); // At P == 1 so k failures should be infinite. + CHECK_EQUAL(quantile(complement(dist, static_cast(0))), std::numeric_limits::infinity()); // Q zero so P == 1. + } else - { // real_concept case, so check it throws rather than returning infinity. - BOOST_CHECK_EQUAL( - quantile( // At P == 1 so k failures should be infinite. - negative_binomial_distribution(static_cast(8), static_cast(0.25)), - static_cast(1)), - boost::math::tools::max_value() ); - - BOOST_CHECK_EQUAL( - quantile(complement( // Q zero 1 so P == 1 < cdf(0) so should be exactly infinity. - negative_binomial_distribution(static_cast(8), static_cast(0.25)), - static_cast(0))), - boost::math::tools::max_value()); + { // real_concept case, which returns max_value rather than infinity. + CHECK_EQUAL(quantile(dist, static_cast(1)), boost::math::tools::max_value()); + CHECK_EQUAL(quantile(complement(dist, static_cast(0))), boost::math::tools::max_value()); } - BOOST_CHECK( // Should work for built-in and real_concept. - quantile(complement( // Q very near to 1 so P nearly 1 < so should be large > 384. - negative_binomial_distribution(static_cast(8), static_cast(0.25)), - static_cast(boost::math::tools::min_value()))) - >= static_cast(384) ); - - BOOST_CHECK_EQUAL( - quantile( // P == 0 < cdf(0) so should be zero. - negative_binomial_distribution(static_cast(8), static_cast(0.25)), - static_cast(0)), - static_cast(0)); + // Q very near to 1 so P nearly 1 so should be large > 384: + CHECK_GE(quantile(complement(dist, boost::math::tools::min_value())), static_cast(384)); + CHECK_EQUAL(quantile(dist, static_cast(0)), static_cast(0)); // P == 0 < cdf(0) so should be zero. // Quantile Complement boundary cases: - - BOOST_CHECK_EQUAL( - quantile(complement( // Q = 1 so P = 0 < cdf(0) so should be exactly zero. - negative_binomial_distribution(static_cast(8), static_cast(0.25)), - static_cast(1))), - static_cast(0) - ); - - BOOST_CHECK_EQUAL( - quantile(complement( // Q very near 1 so P == epsilon < cdf(0) so should be exactly zero. - negative_binomial_distribution(static_cast(8), static_cast(0.25)), - static_cast(1 - boost::math::tools::epsilon()))), - static_cast(0) - ); + CHECK_EQUAL(quantile(complement(dist, static_cast(1))), static_cast(0)); // Q = 1 so P = 0 < cdf(0) + CHECK_EQUAL(quantile(complement(dist, static_cast(1 - boost::math::tools::epsilon()))), static_cast(0)); // P == epsilon < cdf(0) // Check that duff arguments throw domain_error: - BOOST_MATH_CHECK_THROW( - pdf( // Negative successes! - negative_binomial_distribution(static_cast(-1), static_cast(0.25)), - static_cast(0)), std::domain_error - ); - BOOST_MATH_CHECK_THROW( - pdf( // Negative success_fraction! - negative_binomial_distribution(static_cast(8), static_cast(-0.25)), - static_cast(0)), std::domain_error - ); - BOOST_MATH_CHECK_THROW( - pdf( // Success_fraction > 1! - negative_binomial_distribution(static_cast(8), static_cast(1.25)), - static_cast(0)), - std::domain_error - ); - BOOST_MATH_CHECK_THROW( - pdf( // Negative k argument ! - negative_binomial_distribution(static_cast(8), static_cast(0.25)), - static_cast(-1)), - std::domain_error - ); - //BOOST_MATH_CHECK_THROW( - //pdf( // Unlike binomial there is NO limit on k (failures) - //negative_binomial_distribution(static_cast(8), static_cast(0.25)), - //static_cast(9)), std::domain_error - //); - BOOST_MATH_CHECK_THROW( - cdf( // Negative k argument ! - negative_binomial_distribution(static_cast(8), static_cast(0.25)), - static_cast(-1)), - std::domain_error - ); - BOOST_MATH_CHECK_THROW( - cdf( // Negative success_fraction! - negative_binomial_distribution(static_cast(8), static_cast(-0.25)), - static_cast(0)), std::domain_error - ); - BOOST_MATH_CHECK_THROW( - cdf( // Success_fraction > 1! - negative_binomial_distribution(static_cast(8), static_cast(1.25)), - static_cast(0)), std::domain_error - ); - BOOST_MATH_CHECK_THROW( - quantile( // Negative success_fraction! - negative_binomial_distribution(static_cast(8), static_cast(-0.25)), - static_cast(0)), std::domain_error - ); - BOOST_MATH_CHECK_THROW( - quantile( // Success_fraction > 1! - negative_binomial_distribution(static_cast(8), static_cast(1.25)), - static_cast(0)), std::domain_error - ); + CHECK_THROW(pdf(nb(-1, 0.25), static_cast(0)), std::domain_error); // Negative successes! + CHECK_THROW(pdf(nb(8, -0.25), static_cast(0)), std::domain_error); // Negative success_fraction! + CHECK_THROW(pdf(nb(8, 1.25), static_cast(0)), std::domain_error); // Success_fraction > 1! + CHECK_THROW(pdf(nb(8, 0.25), static_cast(-1)), std::domain_error); // Negative k argument! + // Unlike binomial there is NO limit on k (failures). + CHECK_THROW(cdf(nb(8, 0.25), static_cast(-1)), std::domain_error); // Negative k argument! + CHECK_THROW(cdf(nb(8, -0.25), static_cast(0)), std::domain_error); // Negative success_fraction! + CHECK_THROW(cdf(nb(8, 1.25), static_cast(0)), std::domain_error); // Success_fraction > 1! + CHECK_THROW(quantile(nb(8, -0.25), static_cast(0)), std::domain_error); // Negative success_fraction! + CHECK_THROW(quantile(nb(8, 1.25), static_cast(0)), std::domain_error); // Success_fraction > 1! // End of check throwing 'duff' out-of-domain values. +// The table's SC_ uses BOOST_JOIN, from Boost.Config, which a standalone build doesn't have: +#ifndef BOOST_JOIN +#define BOOST_JOIN(X, Y) BOOST_MATH_JOIN(X, Y) +#endif #define T RealType #include "negative_binomial_quantile.ipp" - for(unsigned i = 0; i < negative_binomial_quantile_data.size(); ++i) + for(auto const& row : negative_binomial_quantile_data) { using namespace boost::math::policies; typedef policy > P1; @@ -752,60 +338,53 @@ if(std::numeric_limits::is_specialized) typedef policy > P4; typedef policy > P5; typedef policy > P6; + RealType r = static_cast(row[0]); + RealType sf = static_cast(row[1]); + RealType p = static_cast(row[2]); + RealType lower = static_cast(row[3]); + RealType upper = static_cast(row[4]); RealType tol = boost::math::tools::epsilon() * 700; - if(!boost::is_floating_point::value) + if(!std::is_floating_point::value) tol *= 10; // no lanczos approximation implies less accuracy // - // Check full real value first: + // Check full real value first, to the cdf's tolerance times the quantile's condition number + // (548 for r = 18.3, p = 0.906, where the complement's quantile is 0.0061): // - negative_binomial_distribution p1(negative_binomial_quantile_data[i][0], negative_binomial_quantile_data[i][1]); - RealType x = quantile(p1, negative_binomial_quantile_data[i][2]); - BOOST_CHECK_CLOSE_FRACTION(x, negative_binomial_quantile_data[i][3], tol); - x = quantile(complement(p1, negative_binomial_quantile_data[i][2])); - BOOST_CHECK_CLOSE_FRACTION(x, negative_binomial_quantile_data[i][4], tol); + negative_binomial_distribution p1(r, sf); + CHECK_RELATIVE_ERROR(lower, quantile(p1, p), tol * quantile_condition_number(p1, lower, p, false)); + CHECK_RELATIVE_ERROR(upper, quantile(complement(p1, p)), tol * quantile_condition_number(p1, upper, p, true)); // // Now with round down to integer: // - negative_binomial_distribution p2(negative_binomial_quantile_data[i][0], negative_binomial_quantile_data[i][1]); - x = quantile(p2, negative_binomial_quantile_data[i][2]); - BOOST_CHECK_EQUAL(x, floor(negative_binomial_quantile_data[i][3])); - x = quantile(complement(p2, negative_binomial_quantile_data[i][2])); - BOOST_CHECK_EQUAL(x, floor(negative_binomial_quantile_data[i][4])); + negative_binomial_distribution p2(r, sf); + CHECK_EQUAL(quantile(p2, p), RealType(floor(lower))); + CHECK_EQUAL(quantile(complement(p2, p)), RealType(floor(upper))); // // Now with round up to integer: // - negative_binomial_distribution p3(negative_binomial_quantile_data[i][0], negative_binomial_quantile_data[i][1]); - x = quantile(p3, negative_binomial_quantile_data[i][2]); - BOOST_CHECK_EQUAL(x, ceil(negative_binomial_quantile_data[i][3])); - x = quantile(complement(p3, negative_binomial_quantile_data[i][2])); - BOOST_CHECK_EQUAL(x, ceil(negative_binomial_quantile_data[i][4])); + negative_binomial_distribution p3(r, sf); + CHECK_EQUAL(quantile(p3, p), RealType(ceil(lower))); + CHECK_EQUAL(quantile(complement(p3, p)), RealType(ceil(upper))); // // Now with round to integer "outside": // - negative_binomial_distribution p4(negative_binomial_quantile_data[i][0], negative_binomial_quantile_data[i][1]); - x = quantile(p4, negative_binomial_quantile_data[i][2]); - BOOST_CHECK_EQUAL(x, negative_binomial_quantile_data[i][2] < 0.5f ? floor(negative_binomial_quantile_data[i][3]) : ceil(negative_binomial_quantile_data[i][3])); - x = quantile(complement(p4, negative_binomial_quantile_data[i][2])); - BOOST_CHECK_EQUAL(x, negative_binomial_quantile_data[i][2] < 0.5f ? ceil(negative_binomial_quantile_data[i][4]) : floor(negative_binomial_quantile_data[i][4])); + negative_binomial_distribution p4(r, sf); + CHECK_EQUAL(quantile(p4, p), RealType(p < 0.5f ? floor(lower) : ceil(lower))); + CHECK_EQUAL(quantile(complement(p4, p)), RealType(p < 0.5f ? ceil(upper) : floor(upper))); // // Now with round to integer "inside": // - negative_binomial_distribution p5(negative_binomial_quantile_data[i][0], negative_binomial_quantile_data[i][1]); - x = quantile(p5, negative_binomial_quantile_data[i][2]); - BOOST_CHECK_EQUAL(x, negative_binomial_quantile_data[i][2] < 0.5f ? ceil(negative_binomial_quantile_data[i][3]) : floor(negative_binomial_quantile_data[i][3])); - x = quantile(complement(p5, negative_binomial_quantile_data[i][2])); - BOOST_CHECK_EQUAL(x, negative_binomial_quantile_data[i][2] < 0.5f ? floor(negative_binomial_quantile_data[i][4]) : ceil(negative_binomial_quantile_data[i][4])); + negative_binomial_distribution p5(r, sf); + CHECK_EQUAL(quantile(p5, p), RealType(p < 0.5f ? ceil(lower) : floor(lower))); + CHECK_EQUAL(quantile(complement(p5, p)), RealType(p < 0.5f ? floor(upper) : ceil(upper))); // // Now with round to nearest integer: // - negative_binomial_distribution p6(negative_binomial_quantile_data[i][0], negative_binomial_quantile_data[i][1]); - x = quantile(p6, negative_binomial_quantile_data[i][2]); - BOOST_CHECK_EQUAL(x, floor(negative_binomial_quantile_data[i][3] + 0.5f)); - x = quantile(complement(p6, negative_binomial_quantile_data[i][2])); - BOOST_CHECK_EQUAL(x, floor(negative_binomial_quantile_data[i][4] + 0.5f)); + negative_binomial_distribution p6(r, sf); + CHECK_EQUAL(quantile(p6, p), RealType(floor(lower + 0.5f))); + CHECK_EQUAL(quantile(complement(p6, p)), RealType(floor(upper + 0.5f))); } - - return; +#undef T } // template void test_spots(RealType) // Any floating-point type RealType. template @@ -816,80 +395,53 @@ void test_entropy(RealType) RealType tol = 16 * boost::math::tools::epsilon(); auto check = [&](RealType a, RealType b, RealType h, RealType tolerance) { - BOOST_CHECK_CLOSE_FRACTION(entropy(negative_binomial_distribution(a, b)), h, tolerance); + CHECK_RELATIVE_ERROR(h, entropy(negative_binomial_distribution(a, b)), tolerance); }; check(static_cast(5.0L), static_cast(0.25L), static_cast(3.39359238247812518456836017160705370313434987L), tol); check(static_cast(0.5L), static_cast(0.5L), static_cast(0.943370759269334533806787413559674682998205284L), tol); check(static_cast(1000.5L), static_cast(0.0625L), static_cast(7.61305197413188991442493027236853895416590232L), tol); check(static_cast(1000000.0L), static_cast(0.5L), static_cast(8.67326702746660536495437769998638212050453356L), tol); - // About 1e4 terms; the sum is compensated, which a compiler that reassociates floating-point - // arithmetic (icpx's default) undoes, leaving ~70 eps: - check(static_cast(0.125L), static_cast(0.00390625L), static_cast(3.05357381884703876385483729149700634949203959L), 128 * tol / 16); + // About 1e4 terms, summed with compensation: + check(static_cast(0.125L), static_cast(0.00390625L), static_cast(3.05357381884703876385483729149700634949203959L), tol); check(static_cast(2.5L), static_cast(0.001953125L), static_cast(7.96729439209919695998955366089632316751290691L), tol); check(static_cast(100000000.0L), static_cast(0.9990234375L), static_cast(7.16451918702377261816611719177425949753810124L), tol); - BOOST_CHECK_EQUAL(entropy(negative_binomial_distribution(3, 1)), RealType(0)); + CHECK_EQUAL(entropy(negative_binomial_distribution(3, 1)), RealType(0)); } -BOOST_AUTO_TEST_CASE( test_main ) +int main() { // Check that can generate negative_binomial distribution using the two convenience methods: using namespace boost::math; - negative_binomial mynb1(2., 0.5); // Using typedef - default type is double. - negative_binomial_distribution<> myf2(2., 0.5); // Using default RealType double. - - // Basic sanity-check spot values. - - // Test some simple double only examples. - negative_binomial_distribution my8dist(8., 0.25); - // 8 successes (r), 0.25 success fraction = 35% or 1 in 4 successes. - // Note: double values (matching the distribution definition) avoid the need for any casting. + negative_binomial mynb1(2., 0.5); // Using typedef - default type is double. + negative_binomial_distribution<> myf2(2., 0.5); // Using default RealType double. + static_cast(mynb1); + static_cast(myf2); // Check accessor functions return exact values for double at least. - BOOST_CHECK_EQUAL(my8dist.successes(), static_cast(8)); - BOOST_CHECK_EQUAL(my8dist.success_fraction(), static_cast(1./4.)); + negative_binomial_distribution my8dist(8., 0.25); + CHECK_EQUAL(my8dist.successes(), static_cast(8)); + CHECK_EQUAL(my8dist.success_fraction(), static_cast(1./4.)); // (Parameter value, arbitrarily zero, only communicates the floating point type). #ifdef TEST_FLOAT - test_spots(0.0F); // Test float. + test_spots(0.0F); test_entropy(0.0F); #endif #ifdef TEST_DOUBLE - test_spots(0.0); // Test double. + test_spots(0.0); test_entropy(0.0); #endif #ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS #ifdef TEST_LDOUBLE - test_spots(0.0L); // Test long double. + test_spots(0.0L); test_entropy(0.0L); #endif #ifndef BOOST_MATH_NO_REAL_CONCEPT_TESTS #ifdef TEST_REAL_CONCEPT - test_spots(boost::math::concepts::real_concept(0.)); // Test real concept. - test_entropy(boost::math::concepts::real_concept(0.)); + test_spots(boost::math::concepts::real_concept(0.)); + test_entropy(boost::math::concepts::real_concept(0.)); #endif - #endif -#else - std::cout << "The long double tests have been disabled on this platform " - "either because the long double overloads of the usual math functions are " - "not available at all, or because they are too inaccurate for these tests " - "to pass." << std::endl; #endif - - -} // BOOST_AUTO_TEST_CASE( test_main ) - -/* - -Autorun "i:\boost-06-05-03-1300\libs\math\test\Math_test\debug\test_negative_binomial.exe" -Running 1 test case... -Tolerance = 0.0119209%. -Tolerance 5 eps = 5.96046e-007%. -Tolerance = 2.22045e-011%. -Tolerance 5 eps = 1.11022e-015%. -Tolerance = 2.22045e-011%. -Tolerance 5 eps = 1.11022e-015%. -Tolerance = 2.22045e-011%. -Tolerance 5 eps = 1.11022e-015%. -*** No errors detected - -*/ +#endif + return boost::math::test::report_errors(); +} diff --git a/test/test_out_of_range.hpp b/test/test_out_of_range.hpp index 8b7e723c07..2971ec9708 100644 --- a/test/test_out_of_range.hpp +++ b/test/test_out_of_range.hpp @@ -10,7 +10,12 @@ #include #include -#include +// Included after math_unit_test.hpp, the checks use its CHECK_THROW instead of Boost.Test: +#ifdef BOOST_MATH_TEST_TEST_HPP +# define BOOST_CHECK_THROW(x, y) CHECK_THROW(x, y) +#else +# include +#endif /*` check_out_of_range functions check that bad parameters passed to constructors and functions throw domain_error exceptions. diff --git a/test/test_pow1p.cpp b/test/test_pow1p.cpp index 4e07646122..330f3f4a5d 100644 --- a/test/test_pow1p.cpp +++ b/test/test_pow1p.cpp @@ -68,8 +68,6 @@ void test() CHECK_EQUAL(boost::math::isnan(boost::math::pow1p(boost::math::numeric_limits::quiet_NaN(), boost::math::numeric_limits::infinity())), true); } - // icpx defaults to -fp-model=fast, which breaks the overflow handling and error-free transforms. - #ifndef BOOST_MATH_ENABLE_SYCL // Tiny x, huge finite y: the exponent y*log1p(x) is far out of range, and the // result must overflow or underflow cleanly rather than become inf*0 = NaN. using std::ldexp; @@ -79,7 +77,6 @@ void test() CHECK_EQUAL(boost::math::pow1p(tiny, -huge), T(0)); CHECK_EQUAL(boost::math::pow1p(-tiny, huge), T(0)); CHECK_EQUAL(boost::math::pow1p(-tiny, -huge), boost::math::numeric_limits::infinity()); - #endif // pow(+/-inf, y) CHECK_EQUAL(boost::math::pow1p(boost::math::numeric_limits::infinity(), T(2)), boost::math::numeric_limits::infinity()); @@ -100,7 +97,6 @@ void test() // (1+x) < 0 CHECK_ULP_CLOSE(boost::math::pow1p(T(-3), T(2)), pow(T(-2), T(2)), 10); - #ifndef BOOST_MATH_ENABLE_SYCL // MinGW's long double exp and log lose hundreds of ulp for these huge results, so skip them there. #if defined(__MINGW32__) || defined(__CYGWIN__) const bool check_huge = boost::math::tools::digits() <= 53; @@ -137,7 +133,6 @@ void test() CHECK_ULP_CLOSE(T(1.9141657450800789253871830037475352857803e-261L), boost::math::pow1p(T(-2.1413450492090156e-22), T(2.8034969310889015e+24)), 10); CHECK_ULP_CLOSE(T(1.5108017556765851783235726673009126604833e+287L), boost::math::pow1p(T(6.3377799876470245e-17), T(1.0433536087590337e+19)), 10); } - #endif // x < 0 std::mt19937_64 rng;