Improve precision and speed of the Kolmogorov-Smirnov distribution - #1468
Open
evanmiller wants to merge 2 commits into
Open
evanmiller wants to merge 2 commits into
evanmiller wants to merge 2 commits into
Conversation
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 <noreply@anthropic.com>
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 <noreply@anthropic.com>
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
1 participant
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Measured against a 50-digit reference (
cpp_bin_float_50through the same header), the Kolmogorov-Smirnov distribution had three problems. This PR fixes them without touching the Jacobi theta functions, which the CDF still uses.Quantile
The Newton-Raphson bracket was hard-coded to
[0, 1]inx, so for smalln, where the upper-tail quantile exceeds unity, it either threw (quantile(dist(1), 0.9)) or was clamped: forn = 10the complement quantile of anyq <= 1e-10returned exactly 1. Mid-range it needed about 9 iterations.The iteration now runs on
ln cdfinW = pi^2/(8 x^2 n)forp < 0.6, and onln(1 - cdf)inV = x^2 nabove, where the residual is very nearly linear. The starting value inverts the first two terms of the series in each variable; that is already exact to working precision over most of the range, so a single evaluation usually confirms it.p = 0returns 0 andp = 1raises an overflow error, as in the other distributions.Mode, median, skewness, kurtosis
These are constants of the standardized distribution times a power of
n. The mode was found at run-time by Brent's minimizer, which cannot locate a flat maximum better thansqrt(epsilon); skewness and kurtosis were computed from closed forms that cancel away 40-260 ulp. They are now tabulated to 100 decimal digits (verified two independent ways), with run-time fallbacks for types of higher precision. Amedian()overload is added.PDF
Each series is evaluated from a single exponential, with the remaining terms obtained by multiplication, and the rounding of
x*x*n, ofpi^2/8and of the division is tracked with Dekker products so it is not amplified by the size of the exponent. Just above the smallest normal number the leading term is evaluated as a square ofexp(-W/2), which fixes a loss of three decimal digits to an intermediate underflow (visible as the clipped points in the old ULP graph).Results,
n = 10, doubleA random stress test of 640,000 quantiles (
nfrom 1 to 1e6,plog-uniform over the full range, float and double) had no failures.Tests and docs
test_kolmogorov_smirnov.cppnow covers smalln, both tails down to the smallest normal number (by bracketing, which stays meaningful where one ulp ofxis many ulps ofp), reference values computed at 120 digits, the moment constants against their closed forms, and the high-precision fallbacks with a 110-digit type. The implementation table in the docs is updated and the PDF ULP graph regenerated; its vertical axis is now ±5 ulp.Not addressed: the complement CDF keeps its 10-120 ulp error in the far upper tail, because the rounding of
tauis amplified inside the theta function. That is only fixable inside theta, or by evaluating thez = 0series here directly.🤖 Generated with Claude Code