in commons-statistics-inference/src/main/java/org/apache/commons/statistics/inference/KolmogorovSmirnovDistribution.java [533:635]
static double pelzGood(double x, int n) {
// Change the variable to z since approximation is for the distribution evaluated at d / sqrt(n)
final double z2 = x * x * n;
double lne = -PI2 / (8 * z2);
// Final result is ~ (1 - K0) ~ 1 - sqrt(2pi/z) exp(-pi^2 / (8 z^2))
// Do not compute when the exp value is far below eps.
if (lne < LOG_PG_MIN) {
// z ~ sqrt(-pi^2/(8*min)) ~ 0.1491
return 1;
}
// Note that summing K1, ..., K3 over all k with factor
// (k + 1/2) is equivalent to summing over all k with
// 2 (k - 1/2) / 2 == (2k - 1) / 2
// This is the form for K0.
// Compute all together over odd integers and divide factors
// of (k + 1/2)^b by 2^b.
double k0 = 0;
double k1 = 0;
double k2 = 0;
double k3 = 0;
final double rootN = Math.sqrt(n);
final double z = x * rootN;
final double z3 = z * z2;
final double z4 = z2 * z2;
final double z6 = Math.pow(z2, 3);
final double z7 = Math.pow(z2, 3.5);
final double z8 = Math.pow(z2, 4);
final double z10 = Math.pow(z2, 5);
final double a1 = PI2 / 4;
final double a2 = 6 * z6 + 2 * z4;
final double b2 = (PI2 * (2 * z4 - 5 * z2)) / 4;
final double c2 = (PI4 * (1 - 2 * z2)) / 16;
final double a3 = (PI6 * (5 - 30 * z2)) / 64;
final double b3 = (PI4 * (-60 * z2 + 212 * z4)) / 16;
final double c3 = (PI2 * (135 * z4 - 96 * z6)) / 4;
final double d3 = -(30 * z6 + 90 * z8);
// Iterate j=(2k - 1) for k=1, 2, ...
// Terms reduce in size. Stop when:
// exp(-pi^2 / 8z^2) * eps = exp((2k-1)^2 * -pi^2 / 8z^2)
// (2k-1)^2 = 1 - log(eps) * 8z^2 / pi^2
// 0 = k^2 - k + log(eps) * 2z^2 / pi^2
// Solve using quadratic equation and eps = ulp(1.0): 4a ~ -288
final int max = (int) Math.ceil((1 + Math.sqrt(1 - FOUR_A * z2 / PI2)) / 2);
// Sum smallest terms first
for (int k = max; k > 0; k--) {
final int j = 2 * k - 1;
// Create (2k-1)^2; (2k-1)^4; (2k-1)^6
final double j2 = (double) j * j;
final double j4 = Math.pow(j, 4);
final double j6 = Math.pow(j, 6);
// exp(-pi^2 * (2k-1)^2 / 8z^2)
final double e = Math.exp(lne * j2);
k0 += e;
k1 += (a1 * j2 - z2) * e;
k2 += (a2 + b2 * j2 + c2 * j4) * e;
k3 += (a3 * j6 + b3 * j4 + c3 * j2 + d3) * e;
}
k0 *= ROOT_TWO_PI / z;
// Factors are halved as the sum is for k in -inf to +inf
k1 *= ROOT_HALF_PI / (3 * z4);
k2 *= ROOT_HALF_PI / (36 * z7);
k3 *= ROOT_HALF_PI / (3240 * z10);
// Compute additional K2,K3 terms
double k2b = 0;
double k3b = 0;
// -pi^2 / (2z^2)
lne *= 4;
final double a3b = 3 * PI2 * z2;
// Iterate for j=1, 2, ...
// Note: Here max = sqrt(1 - FOUR_A z^2 / (4 pi^2)).
// This is marginally smaller so we reuse the same value.
for (int j = max; j > 0; j--) {
final double j2 = (double) j * j;
final double j4 = Math.pow(j, 4);
// exp(-pi^2 * k^2 / 2z^2)
final double e = Math.exp(lne * j2);
k2b += PI2 * j2 * e;
k3b += (-PI4 * j4 + a3b * j2) * e;
}
// Factors are halved as the sum is for k in -inf to +inf
k2b *= ROOT_HALF_PI / (18 * z3);
k3b *= ROOT_HALF_PI / (108 * z6);
// Series: K0(z) + K1(z)/n^0.5 + K2(z)/n + K3(z)/n^1.5 + O(1/n^2)
k1 /= rootN;
k2 /= n;
k3 /= n * rootN;
k2b /= n;
k3b /= n * rootN;
// Return (1 - CDF) with an extended precision sum in order of descending magnitude
return clipProbability(Sum.of(1, -k0, -k1, -k2, -k3, +k2b, -k3b).getAsDouble());
}