![]() |
Home | Libraries | People | FAQ | More |
#include <boost/math/distributions/kolmogorov_smirnov.hpp>
namespace boost{ namespace math{ template <class RealType = double, class Policy = policies::policy<> > class kolmogorov_smirnov_distribution; typedef kolmogorov_smirnov_distribution<> kolmogorov_smirnov; template <class RealType, class Policy> class kolmogorov_smirnov_distribution { public: typedef RealType value_type; typedef Policy policy_type; // Constructor: kolmogorov_smirnov_distribution(RealType n); // Accessor to parameter: RealType number_of_observations()const; }} // namespaces
The Kolmogorov-Smirnov test in statistics compares two empirical distributions, or an empirical distribution against any theoretical distribution.[1] It makes use of a specific distribution which is informally known in the literature as the Kolmogorv-Smirnov distribution, implemented here.
Formally, if n observations are taken from a theoretical distribution G(x), and if Gn(x) represents the empirical CDF of those n observations, then the test statistic
will be distributed according to a Kolmogorov-Smirnov distribution parameterized by n.
The exact form of a Kolmogorov-Smirnov distribution is the subject of a large, decades-old literature.[2] In the interest of simplicity, Boost implements the first-order, limiting form of this distribution (the same form originally identified by Kolmogorov[3]), namely
Note that while the exact distribution only has support over [0, 1], this limiting form has positive mass above unity, particularly for small n. The following graph illustrations how the distribution changes for different values of n:
kolmogorov_smirnov_distribution(RealType n);
Constructs a Kolmogorov-Smirnov distribution with n observations.
Requires n > 0, otherwise calls domain_error.
RealType number_of_observations()const;
Returns the parameter n from which this object was constructed.
All the usual non-member accessor functions that are generic to all distributions are supported: Cumulative Distribution Function, Probability Density Function, Quantile, Hazard Function, Cumulative Hazard Function, __logcdf, __logpdf, mean, median, mode, variance, standard deviation, skewness, kurtosis, kurtosis_excess, range and support.
The domain of the random variable is [0, +∞].
The CDF of the Kolmogorov-Smirnov distribution is implemented in terms of the fourth Jacobi Theta function; please refer to the accuracy ULP plots for that function.
The PDF is implemented separately: each of its two series is evaluated from a single exponential whose argument is corrected for the rounding of x2n (and, for small x, of π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:
Because PDF values are simply scaled out and up by the square root of n, the 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.
In the following table, n is the number of observations, x is the random variable, π is Archimedes' constant, and ζ(3) is Apéry's constant.
|
Function |
Implementation Notes |
|---|---|
|
cdf |
When 2*x*x*n <= π: jacobi_theta4tau(0, 2*x*x*n/π) When 2*x*x*n > π: 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 |
|
|
Using a manual derivative of the CDF |
|
cdf complement |
When x*x*n == 0: 1 When 2*x*x*n <= π: 1 - jacobi_theta4tau(0, 2*x*x*n/π) When 2*x*x*n > π: -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 = π2/(8x2n), where ln cdf = ln(4/sqrt(π)) + ln(W)/2 - W + ln(1 + e-8W + e-24W + ...) When p >= 0.6: V = x2n, where ln(1 - cdf) = ln(2) - 2V + ln(1 - e-6V + e-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(π/2) * ln(2) / sqrt(n) |
|
variance |
(π2/12 - π/2*ln2(2))/n |
|
skewness |
0.860426137143668255..., a tabulated evaluation of (9/16*sqrt(π/2)*ζ(3) - 3 * mean * variance - mean3) / (variance3/2) with n = 1 |
|
kurtosis |
3 + 0.881618967910523670..., a tabulated evaluation of (7/720*π4 - 4 * mean * skewness * variance3/2 - 6 * mean2 * variance - mean4) / (variance2) 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.
[2] Simard, R. and L'Ecuyer, P. (2011) "Computing the Two-Sided Kolmogorov-Smirnov Distribution". Journal of Statistical Software, vol. 39, no. 11.
[3] Kolmogorov A (1933). "Sulla determinazione empirica di una legge di distribuzione". G. Ist. Ital. Attuari. 4: 83–91.