From 15e45e0cd256298f0744ba77c831b085232c7adf Mon Sep 17 00:00:00 2001 From: Evan Miller Date: Thu, 17 Sep 2026 06:22:29 -0400 Subject: [PATCH 1/2] Improve precision and speed of the Kolmogorov-Smirnov distribution Quantile: the Newton-Raphson bracket was hard-coded to [0, 1] in x, so for small n the upper-tail quantile (which exceeds unity) either threw or was clamped: quantile(dist(1), 0.9) threw, and for n = 10 the complement quantile of any q <= 1e-10 returned exactly 1. The iteration now runs on the logarithm of the CDF in W = pi^2/(8*x*x*n) for p < 0.6, and on the logarithm of the complement in V = x*x*n above, where the residual is very nearly linear. The starting value inverts the first two terms of the series in each variable, which is already exact to working precision over most of the range, so a single evaluation usually confirms it. p = 0 returns 0 and p = 1 raises an overflow error, as elsewhere. Worst error against a 50-digit reference goes from ~1e6 ulp (1e13 for the complement) to 3 ulp, and the average time roughly halves. Mode, median, skewness and kurtosis are constants of the standardized distribution times a power of n. They are now tabulated to 100 decimal digits (with run-time fallbacks for types of higher precision) instead of being found by Brent's minimizer, which cannot locate a flat maximum better than sqrt(epsilon), or computed from closed forms that cancel away 40-260 ulp. A median() overload is added. PDF: each series is evaluated from a single exponential, the remaining terms coming from multiplication, and the rounding of x*x*n, of pi^2/8 and of the division is tracked with Dekker products so that it is not amplified by the size of the exponent. Mean error drops from 11 to 1.1 ulp and the worst from ~700 to 6 ulp. Just above the smallest normal number the leading term is evaluated as a square of exp(-W/2), which fixes a loss of three decimal digits to an intermediate underflow. The test now covers small n, both tails down to the smallest normal number, reference values computed at 120 digits, the moment constants against their closed forms, and the high-precision fallbacks. The PDF ULP graph is regenerated; the vertical axis is now +-5 ulp. Co-Authored-By: Claude Fable 5.1 --- doc/distributions/kolmogorov_smirnov.qbk | 40 +- doc/graphs/kolmogorov_smirnov_pdf_ulp.svg | 5022 +++++++++-------- .../math/distributions/kolmogorov_smirnov.hpp | 513 +- test/test_kolmogorov_smirnov.cpp | 137 +- 4 files changed, 3056 insertions(+), 2656 deletions(-) diff --git a/doc/distributions/kolmogorov_smirnov.qbk b/doc/distributions/kolmogorov_smirnov.qbk index a748ece697..13f18b98d8 100644 --- a/doc/distributions/kolmogorov_smirnov.qbk +++ b/doc/distributions/kolmogorov_smirnov.qbk @@ -85,8 +85,11 @@ The CDF of the Kolmogorov-Smirnov distribution is implemented in terms of the fourth [link math_toolkit.jacobi_theta.jacobi_theta_overview Jacobi Theta function]; please refer to the accuracy ULP plots for that function. -The PDF is implemented separately, and the following ULP plot illustrates its -accuracy: +The PDF is implemented separately: each of its two series is evaluated from a +single exponential whose argument is corrected for the rounding of x[super 2]n +(and, for small x, of [pi][super 2]/8 and of the division), so that its error +does not grow with the size of the exponent. The following ULP plot illustrates +its accuracy: [graph kolmogorov_smirnov_pdf_ulp] @@ -95,6 +98,12 @@ above plot is representative for all values of /n/. Note that for present purposes, "accuracy" refers to deviations from the limiting approximation, rather than deviations from the exact distribution. +The quantile is accurate to a few ulps of /x/ over the whole range of /p/, +including small /n/ (where it exceeds unity) and probabilities down to the +smallest normalized number. The mode, median, skewness and kurtosis are +tabulated constants of the standardized distribution, accurate to 100 decimal +digits, and so are exact to within the scaling by /n/. + [h4 Implementation] In the following table, /n/ is the number of observations, /x/ is the random variable, @@ -110,15 +119,32 @@ When x*x*n == 0: 1 When 2*x*x*n <= [pi]: 1 - __jacobi_theta4tau(0, 2*x*x*n/[pi]) When 2*x*x*n > [pi]: -__jacobi_theta4m1tau(0, 2*x*x*n/[pi])]] -[[quantile][Using a Newton-Raphson iteration]] -[[quantile from the complement][Using a Newton-Raphson iteration]] -[[mode][Using a run-time PDF maximizer]] +[[quantile][ +Using a Newton-Raphson iteration on the logarithm of the CDF, or of its +complement, in a variable in which it is nearly linear: + +When p < 0.6: W = [pi][super 2]/(8x[super 2]n), where +ln cdf = ln(4/sqrt([pi])) + ln(W)/2 - W + ln(1 + e[super -8W] + e[super -24W] + ...) + +When p >= 0.6: V = x[super 2]n, where +ln(1 - cdf) = ln(2) - 2V + ln(1 - e[super -6V] + e[super -16V] - ...) + +The starting value inverts the first two terms of the expansion, which is +already exact to working precision in the tails.]] +[[quantile from the complement][As the quantile, with p = 1 - q]] +[[mode][0.735467907916571982... / sqrt(n), a tabulated root of the derivative of the PDF]] +[[median][0.827573555189907690... / sqrt(n), a tabulated root of cdf = 1/2]] [[mean][sqrt([pi]/2) * ln(2) / sqrt(n)]] [[variance][([pi][super 2]/12 - [pi]/2*ln[super 2](2))/n]] -[[skewness][(9/16*sqrt([pi]/2)*[zeta](3)/n[super 3/2] - 3 * mean * variance - mean[super 2] * variance) / (variance[super 3/2])]] -[[kurtosis][(7/720*[pi][super 4]/n[super 2] - 4 * mean * skewness * variance[super 3/2] - 6 * mean[super 2] * variance - mean[super 4]) / (variance[super 2])]] +[[skewness][0.860426137143668255..., a tabulated evaluation of (9/16*sqrt([pi]/2)*[zeta](3) - 3 * mean * variance - mean[super 3]) / (variance[super 3/2]) with n = 1]] +[[kurtosis][3 + 0.881618967910523670..., a tabulated evaluation of (7/720*[pi][super 4] - 4 * mean * skewness * variance[super 3/2] - 6 * mean[super 2] * variance - mean[super 4]) / (variance[super 2]) with n = 1]] ] +The tabulated constants are accurate to 100 decimal digits; for types with more +precision than that, the mode is refined by a Newton-Raphson iteration, the +median is computed from the quantile, and the skewness and kurtosis are +computed from the formulae above. + [endsect] [/section:kolmogorov_smirnov_dist Kolmogorov-Smirnov] [/ diff --git a/doc/graphs/kolmogorov_smirnov_pdf_ulp.svg b/doc/graphs/kolmogorov_smirnov_pdf_ulp.svg index f1adef5dc2..3af522200c 100644 --- a/doc/graphs/kolmogorov_smirnov_pdf_ulp.svg +++ b/doc/graphs/kolmogorov_smirnov_pdf_ulp.svg @@ -6,23 +6,23 @@ svg { background-color:black; } Kolmogorov-Smirnov PDF (N=10) ULP plot at float precision - + --75 +-5.582 --50 +-4.109 --25 +-2.637 -0 +-1.164 -25 +0.3087 -50 +1.781 -75 +3.254 -100 +4.727 0.1 @@ -43,2509 +43,2511 @@ svg { background-color:black; } 0.9 1 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/include/boost/math/distributions/kolmogorov_smirnov.hpp b/include/boost/math/distributions/kolmogorov_smirnov.hpp index 66d4997ce8..047460ba17 100644 --- a/include/boost/math/distributions/kolmogorov_smirnov.hpp +++ b/include/boost/math/distributions/kolmogorov_smirnov.hpp @@ -79,15 +79,26 @@ // The PDF is a hand-coded derivative of that function. Actually, there are two // (independent) derivatives, as separate code paths are used for "small x" // (2*x*x*n < pi) and "large x", mirroring the separate code paths in the -// Jacobi Theta implementation to achieve fast convergence. Quantiles are -// computed using a Newton-Raphson iteration from an initial guess that I -// arrived at by trial and error. +// Jacobi Theta implementation to achieve fast convergence. Each path evaluates +// a single exponential and obtains the remaining terms of the series by +// multiplication; the rounding error in the exponent x*x*n is tracked +// explicitly so that it is not amplified by the size of the exponent. +// +// Quantiles are computed by a Newton-Raphson iteration, but not directly in x: +// the lower tail is solved in terms of W = pi^2/(8*x*x*n) and the upper tail +// in terms of V = x*x*n, using the logarithm of the CDF (or its complement) as +// the residual. In these variables the residual is very nearly linear, so the +// iteration converges from anywhere, and the starting point (an inversion of +// the leading terms of the series) is already accurate to a fraction of a +// percent, and to working precision in the tails. // // The mean and variance are implemented using simple closed-form expressions. // Skewness and kurtosis use slightly more complicated closed-form expressions -// that involve the zeta function. The mode is calculated at run-time by -// maximizing the PDF. If you have an analytical solution for the mode, feel -// free to plop it in. +// that involve the zeta function; since these (like the mode and the median) +// do not depend on n except through a factor of sqrt(n), they are stored as +// constants of the standardized (n=1) distribution accurate to 100 decimal +// digits, and only computed at run-time for types with more precision than +// that. // // The CDF and PDF could almost certainly be re-implemented and sped up using a // polynomial or rational approximation, since the only meaningful argument is @@ -111,77 +122,200 @@ #include #include #include +#include +#include #include #include // Newton-Raphson -#include // For the mode namespace boost { namespace math { namespace detail { + +// Constants of the standardized (n=1) distribution, accurate to 100 decimal +// digits. Each accessor scales by the appropriate power of n. Types with more +// precision than the constant (as requested by the policy) fall back to a +// run-time computation. +template +inline bool kolmogorov_smirnov_use_constants(const Policy&) { + return policies::digits() <= 330; +} + +template +inline const RealType& kolmogorov_smirnov_mode_constant() { + static const RealType value = BOOST_MATH_BIG_CONSTANT(RealType, 1000, 0.7354679079165719820624448513051825390913503143340225122926732679283442751250018638342479353073554823); + return value; +} + +template +inline const RealType& kolmogorov_smirnov_median_constant() { + static const RealType value = BOOST_MATH_BIG_CONSTANT(RealType, 1000, 0.8275735551899076901138270828889768075843727832232029452002281054669822747682290224084911530683945353); + return value; +} + +template +inline const RealType& kolmogorov_smirnov_skewness_constant() { + static const RealType value = BOOST_MATH_BIG_CONSTANT(RealType, 1000, 0.8604261371436682558667183685173007452393495203060929924698665373211030237490355834020123201092217303); + return value; +} + +template +inline const RealType& kolmogorov_smirnov_kurtosis_excess_constant() { + static const RealType value = BOOST_MATH_BIG_CONSTANT(RealType, 1000, 0.8816189679105236704015538304328295408355055629289034536289160037385510119639244078254098968815665536); + return value; +} + +// Splits t = hi + lo so that hi carries only the leading half of the +// significand: products of two such hi parts are then exact (Veltkamp's +// splitting, with splitter = 2^ceil(digits/2) + 1). +template +inline RealType kolmogorov_smirnov_splitter() { + BOOST_MATH_STD_USING + return ldexp(RealType(1), (tools::digits() + 1) / 2) + 1; +} + +template +inline void kolmogorov_smirnov_split(RealType t, RealType splitter, RealType& hi, RealType& lo) { + RealType g = splitter * t; + hi = g - (g - t); + lo = t - hi; +} + +// Returns the rounding error of the product ab = a*b, i.e. the exact product +// is ab + result (Dekker's algorithm, which works for any floating-point +// type, unlike fma). +template +inline RealType kolmogorov_smirnov_product_error(RealType a, RealType b, RealType ab, RealType splitter) { + RealType ah, al, bh, bl; + kolmogorov_smirnov_split(a, splitter, ah, al); + kolmogorov_smirnov_split(b, splitter, bh, bl); + return (((ah * bh - ab) + ah * bl) + al * bh) + al * bl; +} + +// Given u = x*x*n as computed in floating point, returns the rounding error +// in u, so that exp(-c*u) can be evaluated as exp(-c*u) * (1 - c*err) +// without the rounding error of u being amplified by the size of the +// exponent c*u. +template +inline RealType kolmogorov_smirnov_x2n_error(RealType x, RealType n, RealType u, RealType splitter) { + RealType x2 = x * x; + RealType e1 = kolmogorov_smirnov_product_error(x, x, x2, splitter); + RealType e2 = kolmogorov_smirnov_product_error(x2, n, u, splitter); + return e2 + e1 * n; // x*x*n == u + err +} + +// Returns pi^2/8 as hi + lo, where hi is exactly representable and lo is +// its complement to about 120 bits, so that W = (pi^2/8)/u can be computed +// without the rounding error of the constant being amplified by W. Types +// with more than 116 bits of precision simply use their own rounding of the +// constant (lo = 0). template -inline RealType kolmogorov_smirnov_quantile_guess(RealType p) { - // Choose a starting point for the Newton-Raphson iteration - if (p > 0.9) - return RealType(1.8) - 5 * (1 - p); - if (p < 0.3) - return p + RealType(0.45); - return p + RealType(0.3); +inline RealType kolmogorov_smirnov_pi_sqr_div_eight(RealType& lo) { + BOOST_MATH_STD_USING + if (tools::digits() > 116) { + lo = 0; + return constants::pi_sqr() / 8; + } + // Five pieces of 24 significant bits each; every one is exact in any + // binary floating-point type with at least 24 bits. + lo = ldexp(RealType(5108270), -47) + ldexp(RealType(15913558), -71) + ldexp(RealType(14839001), -95) + ldexp(RealType(8424474), -119); + return ldexp(RealType(10349030), -23); } -// d/dk (theta2(0, 1/(2*k*k/M_PI))/sqrt(2*k*k*M_PI)) +// d/dx (theta2(0, pi/(2*x*x*n))/sqrt(2*x*x*n)) - valid for all x but converges +// quickly when 2*x*x*n < pi. With u = x*x*n and W = pi^2/(8u): +// +// pdf = sqrt(2*pi*n)/u^2 * SUM_{i>=0} exp(-(2i+1)^2 W) * ((i+1/2)^2 pi^2 - u) template RealType kolmogorov_smirnov_pdf_small_x(RealType x, RealType n, const Policy&) { BOOST_MATH_STD_USING - RealType value = RealType(0), delta = RealType(0), last_delta = RealType(0); RealType eps = policies::get_epsilon(); - int i = 0; RealType pi2 = constants::pi_sqr(); - RealType x2n = x*x*n; - if (x2n*x2n == 0.0) { + RealType u = x * x * n; + if (u == 0) return static_cast(0); + RealType c_lo; + RealType c_hi = kolmogorov_smirnov_pi_sqr_div_eight(c_lo); + RealType W = (c_hi + c_lo) / u; + RealType r = exp(-W); + RealType half = 0; + if (r < tools::min_value()) { + // exp(-W) is subnormal or zero although the result, which is larger + // by a factor of order 1/u^2, may not be: only the leading term + // contributes, and it is evaluated below as a square of exp(-W/2) + // to avoid the intermediate underflow. + half = exp(-W / 2); + if (half == 0) + return static_cast(0); } - while (true) { - delta = exp(-RealType(i+0.5)*RealType(i+0.5)*pi2/(2*x2n)) * (RealType(i+0.5)*RealType(i+0.5)*pi2 - x2n); - - if (delta == 0.0) - break; - - if (last_delta != 0.0 && fabs(delta/last_delta) < eps) - break; + // The rounding errors in the constant, in u and in the division, to + // second order, so that exp(-W) can be corrected by a factor (1 - delta): + // W_exact = W + delta, delta = ((c_hi + c_lo) - W*u - W*err) / u, u_exact = u + err + RealType splitter = kolmogorov_smirnov_splitter(); + RealType err = kolmogorov_smirnov_x2n_error(x, n, u, splitter); + RealType Wu = W * u; + RealType delta = ((((c_hi - Wu) - kolmogorov_smirnov_product_error(W, u, Wu, splitter)) + c_lo) - W * err) / u; + if (half != 0) { + RealType leading = half * (1 - delta / 2) * sqrt(constants::root_two_pi() * sqrt(n) * (pi2 / 4 - u)) / u; + return leading * leading; + } + r *= (1 - delta); // exp(-(W + delta)) + RealType r8 = r * r; + r8 *= r8; + r8 *= r8; - value += delta + delta; - last_delta = delta; + // pw = r^((2i+1)^2), mult = r^(8(i+1)): pw_{i+1} = pw_i * mult_i + RealType sum = 0, term, pw = r, mult = r8; + int i = 0; + do { + term = pw * (RealType(i + 0.5) * RealType(i + 0.5) * pi2 - u); + sum += term; + pw *= mult; + mult *= r8; i++; - } + } while (term > eps * sum); - return value * sqrt(n) * constants::root_half_pi() / (x2n*x2n); + return sum * constants::root_two_pi() * sqrt(n) / (u * u); } -// d/dx (theta4(0, 2*x*x*n/M_PI)) +// d/dx (theta4(0, 2*x*x*n/pi)) - valid for all x but converges quickly when +// 2*x*x*n > pi. With u = x*x*n: +// +// pdf = 8*x*n * SUM_{i>=1} (-1)^(i+1) i^2 exp(-2 i^2 u) template inline RealType kolmogorov_smirnov_pdf_large_x(RealType x, RealType n, const Policy&) { BOOST_MATH_STD_USING - RealType value = RealType(0), delta = RealType(0), last_delta = RealType(0); RealType eps = policies::get_epsilon(); - int i = 1; - while (true) { - delta = 8*x*i*i*exp(-2*i*i*x*x*n); - - if (delta == 0.0) - break; - - if (last_delta != 0.0 && fabs(delta / last_delta) < eps) - break; - - if (i%2 == 0) - delta = -delta; + RealType u = x * x * n; + RealType r = exp(-2 * u); + if (r < tools::min_value()) { + // As above: evaluate the leading term as a square so that it does + // not lose precision to an intermediate underflow. + RealType half = exp(-u); + if (half == 0) + return static_cast(0); + RealType err = kolmogorov_smirnov_x2n_error(x, n, u, kolmogorov_smirnov_splitter()); + RealType leading = half * (1 - err) * sqrt(8 * x * n); + return leading * leading; + } + RealType err = kolmogorov_smirnov_x2n_error(x, n, u, kolmogorov_smirnov_splitter()); + r *= (1 - 2 * err); // exp(-2 (u + err)) + RealType r2 = r * r; - value += delta; - last_delta = delta; + // pw = r^(i^2), odd = r^(2i+1): pw_{i+1} = pw_i * odd_i + RealType sum = 0, term, pw = r, odd = r2 * r; + int i = 1; + do { + term = RealType(i) * RealType(i) * pw; + if (i % 2 == 0) + sum -= term; + else + sum += term; + pw *= odd; + odd *= r2; i++; - } + } while (term > eps * sum); - return value * n; + return 8 * x * n * sum; } } // detail @@ -219,56 +353,161 @@ kolmogorov_smirnov_distribution(RealType)->kolmogorov_smirnov_distribution +inline RealType kolmogorov_smirnov_lower_guess(RealType p) { + BOOST_MATH_STD_USING + RealType L = 2 * constants::ln_two() - constants::log_pi() / 2 - log(p); + RealType W = (L > 1) ? L + log(L) / 2 : RealType(1); + for (int i = 0; i < 3; i++) { + RealType e = exp(-8 * W); + RealType h = W - log(W) / 2 - boost::math::log1p(e) - L; + RealType dh = 1 - 1 / (2 * W) + 8 * e / (1 + e); + W -= h / dh; + } + return W; +} + +// Solves 2 * (t - t^4 + t^9 - t^16) = q for t = exp(-2V), a polynomial +// equation that converges in a few steps from t = q/2 when q <= 0.4, and +// returns the starting value V = -ln(t)/2. The iteration is carried out on +// the ratio r = t / (q/2) so that q may be arbitrarily small. +template +inline RealType kolmogorov_smirnov_upper_guess(RealType q) { + BOOST_MATH_STD_USING + RealType t0 = q / 2; + RealType a = t0 * t0 * t0; + RealType b = a * a * t0 * t0; + RealType c = b * a * a * t0; + RealType r = 1; + for (int i = 0; i < 3; i++) { + RealType r3 = r * r * r; + RealType r8 = r3 * r3 * r * r; + RealType r15 = r8 * r3 * r3 * r; + RealType phi = (r - 1) - a * r3 * r + b * r8 * r - c * r15 * r; + RealType dphi = 1 - 4 * a * r3 + 9 * b * r8 - 16 * c * r15; + r -= phi / dphi; + } + return -(log(q) - constants::ln_two() + boost::math::log1p(r - 1)) / 2; +} + template -struct kolmogorov_smirnov_quantile_functor +struct kolmogorov_smirnov_lower_quantile_functor { - kolmogorov_smirnov_quantile_functor(const boost::math::kolmogorov_smirnov_distribution dist, RealType const& p) - : distribution(dist), prob(p) - { - } - - boost::math::tuple operator()(RealType const& x) - { - RealType fx = cdf(distribution, x) - prob; // Difference cdf - value - to minimize. - RealType dx = pdf(distribution, x); // pdf is 1st derivative. - // return both function evaluation difference f(x) and 1st derivative f'(x). - return boost::math::make_tuple(fx, dx); - } + kolmogorov_smirnov_lower_quantile_functor(const kolmogorov_smirnov_distribution& dist, RealType const& p) + : distribution(dist), log_prob(log(p)), scale(constants::pi_sqr() / (8 * dist.number_of_observations())) + { + } + + boost::math::tuple operator()(RealType const& W) + { + BOOST_MATH_STD_USING + RealType x = sqrt(scale / W); + RealType F = cdf(distribution, x); + RealType f = pdf(distribution, x); + // g(W) = ln F(x(W)) - ln p ; dg/dW = (f/F) * dx/dW, dx/dW = -x/(2W) + return boost::math::make_tuple(log(F) - log_prob, -f * x / (2 * W * F)); + } private: - const boost::math::kolmogorov_smirnov_distribution distribution; - RealType prob; + const kolmogorov_smirnov_distribution& distribution; + RealType log_prob; + RealType scale; }; template -struct kolmogorov_smirnov_complementary_quantile_functor +struct kolmogorov_smirnov_upper_quantile_functor { - kolmogorov_smirnov_complementary_quantile_functor(const boost::math::kolmogorov_smirnov_distribution dist, RealType const& p) - : distribution(dist), prob(p) - { - } - - boost::math::tuple operator()(RealType const& x) - { - RealType fx = cdf(complement(distribution, x)) - prob; // Difference cdf - value - to minimize. - RealType dx = -pdf(distribution, x); // pdf is the negative of the derivative of (1-CDF) - // return both function evaluation difference f(x) and 1st derivative f'(x). - return boost::math::make_tuple(fx, dx); - } + kolmogorov_smirnov_upper_quantile_functor(const kolmogorov_smirnov_distribution& dist, RealType const& q) + : distribution(dist), log_prob(log(q)), n(dist.number_of_observations()) + { + } + + boost::math::tuple operator()(RealType const& V) + { + BOOST_MATH_STD_USING + RealType x = sqrt(V / n); + RealType Q = cdf(complement(distribution, x)); + RealType f = pdf(distribution, x); + // g(V) = ln Q(x(V)) - ln q ; dg/dV = -(f/Q) * dx/dV, dx/dV = x/(2V) + return boost::math::make_tuple(log(Q) - log_prob, -f * x / (2 * V * Q)); + } private: - const boost::math::kolmogorov_smirnov_distribution distribution; - RealType prob; + const kolmogorov_smirnov_distribution& distribution; + RealType log_prob; + RealType n; }; +// Common implementation of quantile() and quantile(complement()): p is the +// lower tail probability and q = 1 - p the upper tail probability, whichever +// was supplied by the caller having been computed from the other. template -struct kolmogorov_smirnov_negative_pdf_functor +RealType kolmogorov_smirnov_quantile_imp(const kolmogorov_smirnov_distribution& dist, RealType p, RealType q, const char* function) { - RealType operator()(RealType const& x) { - if (2*x*x < constants::pi()) { - return -kolmogorov_smirnov_pdf_small_x(x, static_cast(1), Policy()); + BOOST_MATH_STD_USING + RealType n = dist.number_of_observations(); + if (p == 0) + return 0; + if (q == 0) + return policies::raise_overflow_error(function, 0, Policy()); + + // Newton-Raphson in W or V is quadratically convergent from a starting + // point that is accurate to at least 1e-5, so the last step overestimates + // the remaining error by many orders of magnitude and a couple of bits + // can be dropped from the convergence criterion. + const int digits = policies::digits() - 2; + const std::uintmax_t max_root_iterations = policies::get_max_root_iterations(); + std::uintmax_t max_iter = max_root_iterations; + RealType result; + if (p < RealType(0.6)) { + RealType W = kolmogorov_smirnov_lower_guess(p); + // Below a thousand times the smallest normalized number, the CDF can + // no longer be evaluated to full precision; the starting value is + // exact to working precision there anyway. + if (p >= 1000 * tools::min_value()) { + W = tools::newton_raphson_iterate( + kolmogorov_smirnov_lower_quantile_functor(dist, p), + W, RealType(0), tools::max_value(), digits, max_iter); + } else { + max_iter = 0; + } + result = sqrt(constants::pi_sqr() / (8 * n * W)); + } else { + RealType V = kolmogorov_smirnov_upper_guess(q); + if (q >= 1000 * tools::min_value()) { + V = tools::newton_raphson_iterate( + kolmogorov_smirnov_upper_quantile_functor(dist, q), + V, RealType(0), tools::max_value(), digits, max_iter); + } else { + max_iter = 0; } - return -kolmogorov_smirnov_pdf_large_x(x, static_cast(1), Policy()); + result = sqrt(V / n); } -}; + if (max_iter >= max_root_iterations) + { + return policies::raise_evaluation_error(function, "Unable to locate solution in a reasonable time:" // LCOV_EXCL_LINE + " either there is no answer to quantile or the answer is infinite. Current best guess is %1%", result, Policy()); // LCOV_EXCL_LINE + } + return result; +} + } // namespace detail BOOST_MATH_EXPORT template @@ -367,7 +606,6 @@ inline RealType cdf(const complemented2_type inline RealType quantile(const kolmogorov_smirnov_distribution& dist, const RealType& p) { - BOOST_MATH_STD_USING static const char* function = "boost::math::quantile(const kolmogorov_smirnov_distribution<%1%>&, %1%)"; // Error check: RealType error_result; @@ -377,49 +615,23 @@ inline RealType quantile(const kolmogorov_smirnov_distribution if(false == detail::check_df(function, n, &error_result, Policy())) return error_result; - RealType k = detail::kolmogorov_smirnov_quantile_guess(p) / sqrt(n); - const int get_digits = policies::digits();// get digits from policy, - std::uintmax_t max_iter = policies::get_max_root_iterations(); // and max iterations. - - RealType result = tools::newton_raphson_iterate(detail::kolmogorov_smirnov_quantile_functor(dist, p), - k, RealType(0), RealType(1), get_digits, max_iter); - if (max_iter >= policies::get_max_root_iterations()) - { - return policies::raise_evaluation_error(function, "Unable to locate solution in a reasonable time:" // LCOV_EXCL_LINE - " either there is no answer to quantile or the answer is infinite. Current best guess is %1%", result, Policy()); // LCOV_EXCL_LINE - } - return result; + return detail::kolmogorov_smirnov_quantile_imp(dist, p, RealType(1 - p), function); } // quantile BOOST_MATH_EXPORT template inline RealType quantile(const complemented2_type, RealType>& c) { - BOOST_MATH_STD_USING - static const char* function = "boost::math::quantile(const kolmogorov_smirnov_distribution<%1%>&, %1%)"; + static const char* function = "boost::math::quantile(const complemented2_type&, %1%>)"; kolmogorov_smirnov_distribution const& dist = c.dist; RealType n = dist.number_of_observations(); // Error check: RealType error_result; - RealType p = c.param; - - if(false == detail::check_probability(function, p, &error_result, Policy())) + RealType q = c.param; + if(false == detail::check_probability(function, q, &error_result, Policy())) return error_result; if(false == detail::check_df(function, n, &error_result, Policy())) return error_result; - RealType k = detail::kolmogorov_smirnov_quantile_guess(RealType(1-p)) / sqrt(n); - - const int get_digits = policies::digits();// get digits from policy, - std::uintmax_t max_iter = policies::get_max_root_iterations(); // and max iterations. - - RealType result = tools::newton_raphson_iterate( - detail::kolmogorov_smirnov_complementary_quantile_functor(dist, p), - k, RealType(0), RealType(1), get_digits, max_iter); - if (max_iter >= policies::get_max_root_iterations()) - { - return policies::raise_evaluation_error(function, "Unable to locate solution in a reasonable time:" // LCOV_EXCL_LINE - " either there is no answer to quantile or the answer is infinite. Current best guess is %1%", result, Policy()); // LCOV_EXCL_LINE - } - return result; + return detail::kolmogorov_smirnov_quantile_imp(dist, RealType(1 - q), q, function); } // quantile (complemented) BOOST_MATH_EXPORT template @@ -432,10 +644,50 @@ inline RealType mode(const kolmogorov_smirnov_distribution& di if(false == detail::check_df(function, n, &error_result, Policy())) return error_result; - std::pair r = boost::math::tools::brent_find_minima( - detail::kolmogorov_smirnov_negative_pdf_functor(), - static_cast(0), static_cast(1), policies::digits()); - return r.first / sqrt(n); + RealType k = detail::kolmogorov_smirnov_mode_constant(); + if (!detail::kolmogorov_smirnov_use_constants(Policy())) { + // Refine by Newton-Raphson on the stationarity condition of the PDF. + // With u = k*k and a_i = (i+1/2)^2 pi^2, the derivative of the small-x + // PDF series with respect to u vanishes when + // + // G(u) = SUM_{i>=0} exp(-a_i/(2u)) (a_i^2/(2u) - 5 a_i/2 + u) = 0 + RealType eps = policies::get_epsilon(); + RealType pi2 = constants::pi_sqr(); + RealType u = k * k; + for (int iter = 0; iter < 100; iter++) { + RealType G = 0, dG = 0, term; + int i = 0; + do { + RealType a = RealType(i + 0.5) * RealType(i + 0.5) * pi2; + RealType e = exp(-a / (2 * u)); + RealType poly = a * a / (2 * u) - 5 * a / 2 + u; + term = e * poly; + G += term; + dG += e * (a / (2 * u * u) * poly - a * a / (2 * u * u) + 1); + i++; + } while (!(fabs(term) <= eps * fabs(G))); + RealType step = G / dG; + u -= step; + if (fabs(step) <= eps * u) + break; + } + k = sqrt(u); + } + return k / sqrt(n); +} + +BOOST_MATH_EXPORT template +inline RealType median(const kolmogorov_smirnov_distribution& dist) +{ + BOOST_MATH_STD_USING + static const char* function = "boost::math::median(const kolmogorov_smirnov_distribution<%1%>&)"; + RealType n = dist.number_of_observations(); + RealType error_result; + if(false == detail::check_df(function, n, &error_result, Policy())) + return error_result; + if (!detail::kolmogorov_smirnov_use_constants(Policy())) + return quantile(dist, RealType(0.5)); + return detail::kolmogorov_smirnov_median_constant() / sqrt(n); } // Mean and variance come directly from @@ -475,6 +727,8 @@ inline RealType skewness(const kolmogorov_smirnov_distribution RealType error_result; if(false == detail::check_df(function, n, &error_result, Policy())) return error_result; + if (detail::kolmogorov_smirnov_use_constants(Policy())) + return detail::kolmogorov_smirnov_skewness_constant(); RealType ex3 = RealType(0.5625) * constants::root_half_pi() * constants::zeta_three() / n / sqrt(n); RealType mean = boost::math::mean(dist); RealType var = boost::math::variance(dist); @@ -482,30 +736,27 @@ inline RealType skewness(const kolmogorov_smirnov_distribution } BOOST_MATH_EXPORT template -inline RealType kurtosis(const kolmogorov_smirnov_distribution& dist) +inline RealType kurtosis_excess(const kolmogorov_smirnov_distribution& dist) { BOOST_MATH_STD_USING - static const char* function = "boost::math::kurtosis(const kolmogorov_smirnov_distribution<%1%>&)"; + static const char* function = "boost::math::kurtosis_excess(const kolmogorov_smirnov_distribution<%1%>&)"; RealType n = dist.number_of_observations(); RealType error_result; if(false == detail::check_df(function, n, &error_result, Policy())) return error_result; + if (detail::kolmogorov_smirnov_use_constants(Policy())) + return detail::kolmogorov_smirnov_kurtosis_excess_constant(); RealType ex4 = 7 * constants::pi_sqr_div_six() * constants::pi_sqr_div_six() / 20 / n / n; RealType mean = boost::math::mean(dist); RealType var = boost::math::variance(dist); RealType skew = boost::math::skewness(dist); - return (ex4 - 4 * mean * skew * var * sqrt(var) - 6 * mean * mean * var - mean * mean * mean * mean) / var / var; + return (ex4 - 4 * mean * skew * var * sqrt(var) - 6 * mean * mean * var - mean * mean * mean * mean) / var / var - 3; } BOOST_MATH_EXPORT template -inline RealType kurtosis_excess(const kolmogorov_smirnov_distribution& dist) +inline RealType kurtosis(const kolmogorov_smirnov_distribution& dist) { - static const char* function = "boost::math::kurtosis_excess(const kolmogorov_smirnov_distribution<%1%>&)"; - RealType n = dist.number_of_observations(); - RealType error_result; - if(false == detail::check_df(function, n, &error_result, Policy())) - return error_result; - return kurtosis(dist) - 3; + return kurtosis_excess(dist) + 3; } }} #endif diff --git a/test/test_kolmogorov_smirnov.cpp b/test/test_kolmogorov_smirnov.cpp index 58dad8cd5d..894c1269c1 100644 --- a/test/test_kolmogorov_smirnov.cpp +++ b/test/test_kolmogorov_smirnov.cpp @@ -12,6 +12,8 @@ #include // for BOOST_CHECK_CLOSE #include #include +#include +#include template // Any floating-point type RealType. void test_spots(RealType) @@ -35,13 +37,14 @@ void test_spots(RealType) RealType x = 1.0 * (i+1) / 1001; RealType p = cdf(dist, x); RealType p1 = cdf(complement(dist, x)); - RealType x1; - if (p < 0.5) - x1 = quantile(dist, p); - else - x1 = quantile(complement(dist, p1)); - if (p > tol && p1 > tol) // skip the extreme tails + if (p > tol && p1 > tol) { // skip the extreme tails + RealType x1; + if (p < 0.5) + x1 = quantile(dist, p); + else + x1 = quantile(complement(dist, p1)); BOOST_CHECK_CLOSE_FRACTION(x, x1, tol); + } } } @@ -87,17 +90,135 @@ void test_spots(RealType) BOOST_CHECK_CLOSE_FRACTION(kurt, kurtosis_excess(dist) + 3, eps); } +// Checks that x is the quantile of p to within a few ulps, by bracketing: +// cdf(x * (1 - margin)) <= p <= cdf(x * (1 + margin)). Unlike a round-trip +// comparison of cdf(quantile(p)) with p, this remains meaningful deep in the +// tails, where the CDF is so steep that one ulp of x corresponds to many +// ulps of p. +template +void check_quantile(const boost::math::kolmogorov_smirnov_distribution& dist, RealType p, RealType x, bool complemented) +{ + using namespace boost::math; + RealType margin = 8 * tools::epsilon(); + RealType lo = complemented ? cdf(complement(dist, RealType(x * (1 + margin)))) : cdf(dist, RealType(x * (1 - margin))); + RealType hi = complemented ? cdf(complement(dist, RealType(x * (1 - margin)))) : cdf(dist, RealType(x * (1 + margin))); + BOOST_TEST_CHECK(lo <= p); + BOOST_TEST_CHECK(p <= hi); +} + +template +void test_quantile_extremes(RealType) +{ + using namespace boost::math; + RealType eps = tools::epsilon(); + + // Small n, where the quantile exceeds unity in the upper tail + for (RealType n : { RealType(1), RealType(2), RealType(3), RealType(1000000) }) { + kolmogorov_smirnov_distribution dist(n); + for (int i = 1; i < 100; i++) { + RealType p = RealType(i) / 100; + check_quantile(dist, p, quantile(dist, p), false); + check_quantile(dist, p, quantile(complement(dist, p)), true); + } + } + + kolmogorov_smirnov_distribution dist(10); + + // Deep tails, down to where the CDF itself is no longer evaluated to + // full precision because its terms become denormal + for (RealType p = RealType(0.1); p > 10000 * tools::min_value(); p /= 1000) { + check_quantile(dist, p, quantile(dist, p), false); + check_quantile(dist, p, quantile(complement(dist, p)), true); + } + + // Limits of the probability range + BOOST_CHECK_EQUAL(quantile(dist, RealType(0)), RealType(0)); + BOOST_CHECK_EQUAL(quantile(complement(dist, RealType(1))), RealType(0)); + BOOST_CHECK_THROW(quantile(dist, RealType(1)), std::overflow_error); + BOOST_CHECK_THROW(quantile(complement(dist, RealType(0))), std::overflow_error); + + // The median is the quantile of one half + BOOST_CHECK_CLOSE_FRACTION(cdf(dist, median(dist)), RealType(0.5), 10 * eps); + BOOST_CHECK_CLOSE_FRACTION(median(dist), quantile(dist, RealType(0.5)), 10 * eps); +} + +// Reference values computed with the same series at 120 decimal digits. +template +void test_reference_values(RealType) +{ + using namespace boost::math; + RealType eps = tools::epsilon(); + static const struct { double n, x; const char *pdf, *cdf, *ccdf; } data[] = { + { 10, 0.125, "0.279400023454107971738809672357", "0.00236117478750956819598968016445", "0.997638825212490431804010319836" }, + { 10, 0.25, "5.19340083129894262893339472019", "0.440440289804736223158552770951", "0.559559710195263776841447229049" }, + { 10, 0.375, "1.80007916388412209267020060796", "0.879916678784370975806015389694", "0.120083321215629024193984610306" }, + { 10, 0.5, "0.269517550178839104001256352353", "0.986524110124136310626593187475", "0.0134758898758636893734068125254" }, + { 10, 0.75, "0.000780437859244050388721178490249", "0.999973985404691864815291512551", "2.60145953081351847084874493052e-05" }, + { 10, 1.0, "1.64892289795084626237275224637e-07", "0.999999995877692755122884344068", "4.12230724487711565593188072421e-09" }, + { 1, 0.25, "4.12855106473622609508308670705e-06", "2.68238100848298275381058297367e-08", "0.999999973176189915170172461894" }, + { 1, 0.5, "0.639582850940456634645459338412", "0.0360547563351249056140861037179", "0.963945243664875094385913896282" }, + { 1, 0.75, "1.6834609513049753816942330657", "0.372832958223738358506334471397", "0.627167041776261641493665528603" }, + { 1, 1.5, "0.133307227419880210037320465613", "0.97778203738347487127945638539", "0.0222179626165251287205436146105" }, + { 1, 2.5, "7.45330634415734044284982458077e-05", "0.999992546693655842658399900267", "7.45330634415734160009973335912e-06" }, + }; + for (const auto& d : data) { + kolmogorov_smirnov_distribution dist(static_cast(d.n)); + RealType x = static_cast(d.x); + BOOST_CHECK_CLOSE_FRACTION(pdf(dist, x), boost::lexical_cast(d.pdf), 20 * eps); + BOOST_CHECK_CLOSE_FRACTION(cdf(dist, x), boost::lexical_cast(d.cdf), 50 * eps); + // The complement is computed by the Jacobi Theta function from a + // rounded argument, whose error is amplified by the exponent 2*x*x*n. + BOOST_CHECK_CLOSE_FRACTION(cdf(complement(dist, x)), boost::lexical_cast(d.ccdf), 200 * eps); + } + + // The moment constants agree with the closed forms they were computed from + kolmogorov_smirnov_distribution dist(10); + RealType n = 10; + RealType mean = boost::math::mean(dist); + RealType var = variance(dist); + RealType ex3 = RealType(0.5625) * constants::root_half_pi() * constants::zeta_three() / n / sqrt(n); + RealType ex4 = 7 * constants::pi_sqr_div_six() * constants::pi_sqr_div_six() / 20 / n / n; + RealType skew = (ex3 - 3 * mean * var - mean * mean * mean) / var / sqrt(var); + RealType kurt = (ex4 - 4 * mean * skew * var * sqrt(var) - 6 * mean * mean * var - mean * mean * mean * mean) / var / var; + BOOST_CHECK_CLOSE_FRACTION(skewness(dist), skew, 500 * eps); + BOOST_CHECK_CLOSE_FRACTION(kurtosis(dist), kurt, 500 * eps); + BOOST_CHECK_CLOSE_FRACTION(kurtosis_excess(dist), kurt - 3, 2000 * eps); +} + +// Types with more than 100 decimal digits compute the mode, median and +// moments at run-time; check that they agree with the tabulated constants. +void test_high_precision() +{ + using namespace boost::math; + typedef boost::multiprecision::number > RealType; + kolmogorov_smirnov_distribution dist(1); + RealType tol = RealType("1e-99"); + BOOST_CHECK_CLOSE_FRACTION(mode(dist), RealType("0.7354679079165719820624448513051825390913503143340225122926732679283442751250018638342479353073554823"), tol); + BOOST_CHECK_CLOSE_FRACTION(median(dist), RealType("0.8275735551899076901138270828889768075843727832232029452002281054669822747682290224084911530683945353"), tol); + BOOST_CHECK_CLOSE_FRACTION(skewness(dist), RealType("0.8604261371436682558667183685173007452393495203060929924698665373211030237490355834020123201092217303"), tol); + BOOST_CHECK_CLOSE_FRACTION(kurtosis_excess(dist), RealType("0.8816189679105236704015538304328295408355055629289034536289160037385510119639244078254098968815665536"), tol); +} + BOOST_AUTO_TEST_CASE( test_main ) { BOOST_MATH_CONTROL_FP; // (Parameter value, arbitrarily zero, only communicates the floating point type). - test_spots(0.0F); // Test float. - test_spots(0.0); // Test double. + test_spots(0.0F); // Test float. + test_spots(0.0); // Test double. + test_quantile_extremes(0.0F); + test_quantile_extremes(0.0); + test_reference_values(0.0F); + test_reference_values(0.0); #ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS test_spots(0.0L); // Test long double. + test_quantile_extremes(0.0L); + test_reference_values(0.0L); #if !defined(BOOST_MATH_NO_REAL_CONCEPT_TESTS) test_spots(boost::math::concepts::real_concept(0.)); // Test real concept. + test_quantile_extremes(boost::math::concepts::real_concept(0.)); + test_reference_values(boost::math::concepts::real_concept(0.)); #endif #endif + test_high_precision(); } From 1d4ba26aa12f62b8dda34fb90144057203a01b7d Mon Sep 17 00:00:00 2001 From: Evan Miller Date: Thu, 17 Sep 2026 11:32:49 -0400 Subject: [PATCH 2/2] Evaluate the upper tail of the Kolmogorov-Smirnov CDF from the nome For 2*x*x*n > pi the CDF and its complement called jacobi_theta4m1tau with tau = 2*x*x*n/pi. The theta function multiplies tau back by pi and by n^2 before exponentiating, so the four roundings in forming tau cost about 4*x*x*n ulps in the result: against a 50-digit reference the complement was off by up to 93 ulps (mean 13) for x*sqrt(n) in [1.25, 6], and no change inside the theta function can recover what the caller has already rounded away. Pass the nome q = exp(-2*x*x*n) instead, with the rounding of x*x*n compensated by the same product-error helper the PDF uses. Below exp(-pi) the theta function raises q to integer powers, which amplifies nothing, and the complement is now within 1.8 ulps (mean 0.5) over the same range. The lower branch (2*x*x*n <= pi) still goes through tau, where it has no equivalent route. The reference-value test checks the complement to 10 epsilon on the upper branch instead of 200. The upper branch costs about 20 ns more per evaluation in double, because the theta function evaluates q^(n^2) with pow() rather than exp(). Co-Authored-By: Claude Fable 5.1 --- doc/distributions/kolmogorov_smirnov.qbk | 7 +++-- .../math/distributions/kolmogorov_smirnov.hpp | 31 +++++++++++++++++-- test/test_kolmogorov_smirnov.cpp | 9 ++++-- 3 files changed, 40 insertions(+), 7 deletions(-) diff --git a/doc/distributions/kolmogorov_smirnov.qbk b/doc/distributions/kolmogorov_smirnov.qbk index 13f18b98d8..58e3c4a7fd 100644 --- a/doc/distributions/kolmogorov_smirnov.qbk +++ b/doc/distributions/kolmogorov_smirnov.qbk @@ -111,14 +111,17 @@ In the following table, /n/ is the number of observations, /x/ is the random var [table [[Function][Implementation Notes]] -[[cdf][Using the relation: cdf = __jacobi_theta4tau(0, 2*x*x/[pi])]] +[[cdf][ +When 2*x*x*n <= [pi]: __jacobi_theta4tau(0, 2*x*x*n/[pi]) + +When 2*x*x*n > [pi]: 1 + __jacobi_theta4m1(0, q) with the nome q = exp(-2*x*x*n), whose exponent is formed with its rounding error compensated so that the theta function, which raises q to integer powers, does not amplify it]] [[pdf][Using a manual derivative of the CDF]] [[cdf complement][ When x*x*n == 0: 1 When 2*x*x*n <= [pi]: 1 - __jacobi_theta4tau(0, 2*x*x*n/[pi]) -When 2*x*x*n > [pi]: -__jacobi_theta4m1tau(0, 2*x*x*n/[pi])]] +When 2*x*x*n > [pi]: -__jacobi_theta4m1(0, q), as for the cdf]] [[quantile][ Using a Newton-Raphson iteration on the logarithm of the CDF, or of its complement, in a variable in which it is nearly linear: diff --git a/include/boost/math/distributions/kolmogorov_smirnov.hpp b/include/boost/math/distributions/kolmogorov_smirnov.hpp index 047460ba17..c90b567499 100644 --- a/include/boost/math/distributions/kolmogorov_smirnov.hpp +++ b/include/boost/math/distributions/kolmogorov_smirnov.hpp @@ -203,6 +203,22 @@ inline RealType kolmogorov_smirnov_x2n_error(RealType x, RealType n, RealType u, return e2 + e1 * n; // x*x*n == u + err } +// The nome q = exp(-2*x*x*n) of the theta function in the CDF, with the +// rounding of x*x*n compensated. Passing q rather than tau = 2*x*x*n/pi +// matters when 2*x*x*n > pi: the theta function then raises q to integer +// powers, so nothing amplifies the rounding, whereas any rounding of tau is +// multiplied back by the size of the exponent. +template +inline RealType kolmogorov_smirnov_nome(RealType x, RealType n) { + BOOST_MATH_STD_USING + RealType u = x * x * n; + RealType q = exp(-2 * u); + if (q == 0) + return q; + RealType err = kolmogorov_smirnov_x2n_error(x, n, u, kolmogorov_smirnov_splitter()); + return q * (1 - 2 * err); // exp(-2 (u + err)) +} + // Returns pi^2/8 as hi + lo, where hi is exactly representable and lo is // its complement to about 120 bits, so that W = (pi^2/8)/u can be computed // without the rounding error of the constant being amplified by W. Types @@ -573,6 +589,13 @@ inline RealType cdf(const kolmogorov_smirnov_distribution& dis if (x*x*n == 0) return 0; + if (2*x*x*n > constants::pi()) { + RealType q = detail::kolmogorov_smirnov_nome(x, n); + if (q == 0) + return 1; + return RealType(1) + jacobi_theta4m1(RealType(0), q, Policy()); + } + return jacobi_theta4tau(RealType(0), 2*x*x*n/constants::pi(), Policy()); } // cdf @@ -597,8 +620,12 @@ inline RealType cdf(const complemented2_type constants::pi()) - return -jacobi_theta4m1tau(RealType(0), 2*x*x*n/constants::pi(), Policy()); + if (2*x*x*n > constants::pi()) { + RealType q = detail::kolmogorov_smirnov_nome(x, n); + if (q == 0) + return 0; + return -jacobi_theta4m1(RealType(0), q, Policy()); + } return RealType(1) - jacobi_theta4tau(RealType(0), 2*x*x*n/constants::pi(), Policy()); } // cdf (complemented) diff --git a/test/test_kolmogorov_smirnov.cpp b/test/test_kolmogorov_smirnov.cpp index 894c1269c1..4b986ce66f 100644 --- a/test/test_kolmogorov_smirnov.cpp +++ b/test/test_kolmogorov_smirnov.cpp @@ -166,9 +166,12 @@ void test_reference_values(RealType) RealType x = static_cast(d.x); BOOST_CHECK_CLOSE_FRACTION(pdf(dist, x), boost::lexical_cast(d.pdf), 20 * eps); BOOST_CHECK_CLOSE_FRACTION(cdf(dist, x), boost::lexical_cast(d.cdf), 50 * eps); - // The complement is computed by the Jacobi Theta function from a - // rounded argument, whose error is amplified by the exponent 2*x*x*n. - BOOST_CHECK_CLOSE_FRACTION(cdf(complement(dist, x)), boost::lexical_cast(d.ccdf), 200 * eps); + // Above 2*x*x*n = pi the complement is computed from the nome with a + // compensated exponent; below it goes through the Jacobi Theta + // function with a rounded tau, whose error is amplified by the + // exponent. + RealType ccdf_tol = (2 * x * x * d.n > constants::pi()) ? 10 * eps : 200 * eps; + BOOST_CHECK_CLOSE_FRACTION(cdf(complement(dist, x)), boost::lexical_cast(d.ccdf), ccdf_tol); } // The moment constants agree with the closed forms they were computed from