Boost C++ Libraries Home Libraries People FAQ More

PrevUpHomeNext

Kolmogorov-Smirnov Distribution

#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:

Member Functions
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.

Non-member Accessors

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, +∞].

Accuracy

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.

Implementation

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

pdf

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.



[3] Kolmogorov A (1933). "Sulla determinazione empirica di una legge di distribuzione". G. Ist. Ital. Attuari. 4: 83–91.


PrevUpHomeNext