diff --git a/doc/distributions/kolmogorov_smirnov.qbk b/doc/distributions/kolmogorov_smirnov.qbk
index a748ece697..58e3c4a7fd 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,
@@ -102,23 +111,43 @@ 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])]]
-[[quantile][Using a Newton-Raphson iteration]]
-[[quantile from the complement][Using a Newton-Raphson iteration]]
-[[mode][Using a run-time PDF maximizer]]
+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:
+
+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..c90b567499 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,216 @@
#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 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 void kolmogorov_smirnov_split(RealType t, RealType splitter, RealType& hi, RealType& lo) {
+ RealType g = splitter * t;
+ hi = g - (g - t);
+ lo = t - hi;
}
-// d/dk (theta2(0, 1/(2*k*k/M_PI))/sqrt(2*k*k*M_PI))
+// 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
+}
+
+// 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
+// with more than 116 bits of precision simply use their own rounding of the
+// constant (lo = 0).
+template
+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/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 +369,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
@@ -334,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
@@ -358,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)
@@ -367,7 +633,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 +642,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 +671,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 +754,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 +763,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..4b986ce66f 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,138 @@ 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);
+ // 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
+ 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();
}