package jhplot.math;
/*
* Class: KolmogorovSmirnovDist
* Environment: Java
* Author: Richard Simard
* Organization: DIRO, Université de Montréal
* Date: 14 december 2010
* Version: 1.0
* Copyright 1 March 2010 by Université de Montréal,
Richard Simard and Pierre L'Ecuyer
=====================================================================
This program is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, version 3 of the License.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with this program. If not, see This class computes both the cumulative probability P[D_n <= x]
* and the complementary cumulative probability P[D_n >= x] of the
* 2-sided 1-sample Kolmogorov-Smirnov distribution.
*
* The Kolmogorov-Smirnov test statistic D_n is defined by
*
* D_n = sup_x |F(x) - S_n(x)|
*
* where n is the sample size, S_n(x) is an empirical distribution function,
* and F(x) is a completely specified theoretical distribution.
*
*
* @author Richard Simard
* @version 1.0
* @since 1 March 2010
*/
public class KolmogorovSmirnovDist
{
private static final double num_Ln2 = 0.69314718055994530941; // ln(2)
private static final double PI2 = Math.PI * Math.PI;
private static final double PI4 = PI2 * PI2;
private static final int NFACT = 20;
private static final int MFACT = 30;
/* For x close to 0 or 1, we use the exact formulae of Ruben-Gambino in all
cases. For n <= NEXACT, we use exact algorithms: the Durbin matrix and
the Pomeranz algorithms. For n > NEXACT, we use asymptotic methods
except for x close to 0 where we still use the method of Durbin
for n <= NKOLMO. For n > NKOLMO, we use asymptotic methods only and
so the precision is less for x close to 0.
We could increase the limit NKOLMO to 10^6 to get better precision
for x close to 0, but at the price of a slower speed. */
private static final int NEXACT = 140;
private static final int NKOLMO = 100000;
// For the Durbin matrix algorithm
private static final double NORM = 1.0e140;
private static final double INORM = 1.0e-140;
private static final int LOGNORM = 140;
//========================================================================
// The factorial n! for 0 <= n <= NFACT
private static final double[] Factorial = {
1,
1,
2,
6,
24,
120,
720,
5040,
40320,
362880,
3628800,
39916800,
479001600.,
6227020800.,
87178291200.,
1307674368000.,
20922789888000.,
355687428096000.,
6402373705728000.,
1.21645100408832e+17,
2.43290200817664e+18
};
//========================================================================
// The natural logarithm of factorial n! for 0 <= n <= MFACT
private static final double[] LnFactorial = {
0.,
0.,
0.6931471805599453,
1.791759469228055,
3.178053830347946,
4.787491742782046,
6.579251212010101,
8.525161361065415,
10.60460290274525,
12.80182748008147,
15.10441257307552,
17.50230784587389,
19.98721449566188,
22.55216385312342,
25.19122118273868,
27.89927138384088,
30.67186010608066,
33.50507345013688,
36.39544520803305,
39.33988418719949,
42.33561646075348,
45.3801388984769,
48.47118135183522,
51.60667556776437,
54.7847293981123,
58.00360522298051,
61.26170176100199,
64.55753862700632,
67.88974313718154,
71.257038967168,
74.65823634883016
};
//------------------------------------------------------------------------
private static double getLogFactorial (int n)
{
// Returns the natural logarithm of factorial n!
if (n <= MFACT) {
return LnFactorial[n];
} else {
double x = (double) (n + 1);
double y = 1.0 / (x * x);
double z = ((-(5.95238095238E-4 * y) + 7.936500793651E-4) * y -
2.7777777777778E-3) * y + 8.3333333333333E-2;
z = ((x - 0.5) * Math.log (x) - x) + 9.1893853320467E-1 + z / x;
return z;
}
}
//========================================================================
private static double KSPlusbarAsymp (int n, double x)
{
/* Compute the probability of the KS+ distribution using an asymptotic
formula */
double t = (6.0 * n * x + 1);
double z = t * t / (18.0 * n);
double v = 1.0 - (2.0 * z * z - 4.0 * z - 1.0) / (18.0 * n);
if (v <= 0.0)
return 0.0;
v = v * Math.exp (-z);
if (v >= 1.0)
return 1.0;
return v;
}
//-------------------------------------------------------------------------
private static double KSPlusbarUpper (int n, double x)
{
/* Compute the probability of the complementary KS+ distribution in
the upper tail using Smirnov's stable formula */
if (n > 200000)
return KSPlusbarAsymp (n, x);
int jmax = (int) (n * (1.0 - x));
// Avoid log(0) for j = jmax and q ~ 1.0
if ((1.0 - x - (double) jmax / n) <= 0.0)
jmax--;
int jdiv;
if (n > 3000)
jdiv = 2;
else
jdiv = 3;
int j = jmax / jdiv + 1;
double LogCom = getLogFactorial (n) - getLogFactorial (j) -
getLogFactorial (n - j);
double LOGJMAX = LogCom;
final double EPSILON = 1.0E-12;
double q;
double term;
double t;
double Sum = 0.0;
while (j <= jmax) {
q = (double) j / n + x;
term = LogCom + (j - 1) * Math.log (q) + (n - j) * Math.log1p (-q);
t = Math.exp (term);
Sum += t;
LogCom += Math.log ((double) (n - j) / (j + 1));
if (t <= Sum * EPSILON)
break;
j++;
}
j = jmax / jdiv;
LogCom = LOGJMAX + Math.log ((double) (j + 1) / (n - j));
while (j > 0) {
q = (double) j / n + x;
term = LogCom + (j - 1) * Math.log (q) + (n - j) * Math.log1p (-q);
t = Math.exp (term);
Sum += t;
LogCom += Math.log ((double) j / (n - j + 1));
if (t <= Sum * EPSILON)
break;
j--;
}
Sum *= x;
// add the term j = 0
Sum += Math.exp (n * Math.log1p (-x));
return Sum;
}
//========================================================================
private static double Pelz (int n, double x)
{
/* Approximating the Lower Tail-Areas of the Kolmogorov-Smirnov One-Sample
Statistic,
Wolfgang Pelz and I. J. Good,
Journal of the Royal Statistical Society, Series B.
Vol. 38, No. 2 (1976), pp. 152-156
*/
final int JMAX = 20;
final double EPS = 1.0e-10;
final double C = 2.506628274631001; // sqrt(2*Pi)
final double C2 = 1.2533141373155001; // sqrt(Pi/2)
final double RACN = Math.sqrt ((double) n);
final double z = RACN * x;
final double z2 = z * z;
final double z4 = z2 * z2;
final double z6 = z4 * z2;
final double w = PI2 / (2.0 * z * z);
double ti, term, tom;
double sum;
int j;
term = 1;
j = 0;
sum = 0;
while (j <= JMAX && term > EPS * sum) {
ti = j + 0.5;
term = Math.exp (-ti * ti * w);
sum += term;
j++;
}
sum *= C / z;
term = 1;
tom = 0;
j = 0;
while (j <= JMAX && Math.abs (term) > EPS * Math.abs (tom)) {
ti = j + 0.5;
term = (PI2 * ti * ti - z2) * Math.exp (-ti * ti * w);
tom += term;
j++;
}
sum += tom * C2 / (RACN * 3.0 * z4);
term = 1;
tom = 0;
j = 0;
while (j <= JMAX && Math.abs (term) > EPS * Math.abs (tom)) {
ti = j + 0.5;
term = 6 * z6 + 2 * z4 + PI2 * (2 * z4 - 5 * z2) * ti * ti +
PI4 * (1 - 2 * z2) * ti * ti * ti * ti;
term *= Math.exp (-ti * ti * w);
tom += term;
j++;
}
sum += tom * C2 / (n * 36.0 * z * z6);
term = 1;
tom = 0;
j = 1;
while (j <= JMAX && term > EPS * tom) {
ti = j;
term = PI2 * ti * ti * Math.exp (-ti * ti * w);
tom += term;
j++;
}
sum -= tom * C2 / (n * 18.0 * z * z2);
term = 1;
tom = 0;
j = 0;
while (j <= JMAX && Math.abs (term) > EPS * Math.abs (tom)) {
ti = j + 0.5;
ti = ti * ti;
term = -30 * z6 - 90 * z6 * z2 + PI2 * (135 * z4 - 96 * z6) * ti +
PI4 * (212 * z4 - 60 * z2) * ti * ti + PI2 * PI4 * ti * ti * ti * (5 -
30 * z2);
term *= Math.exp (-ti * w);
tom += term;
j++;
}
sum += tom * C2 / (RACN * n * 3240.0 * z4 * z6);
term = 1;
tom = 0;
j = 1;
while (j <= JMAX && Math.abs (term) > EPS * Math.abs (tom)) {
ti = j * j;
term = (3 * PI2 * ti * z2 - PI4 * ti * ti) * Math.exp (-ti * w);
tom += term;
j++;
}
sum += tom * C2 / (RACN * n * 108.0 * z6);
return sum;
}
//=========================================================================
private static void CalcFloorCeil (
int n, // sample size
double t, // = nx
double[] A, // A_i
double[] Atflo, // floor (A_i - t)
double[] Atcei // ceiling (A_i + t)
)
{
// Precompute A_i, floors, and ceilings for limits of sums in the
// Pomeranz algorithm
int i;
int ell = (int) t; // floor (t)
double z = t - ell; // t - floor (t)
double w = Math.ceil (t) - t;
if (z > 0.5) {
for (i = 2; i <= 2 * n + 2; i += 2)
Atflo[i] = i / 2 - 2 - ell;
for (i = 1; i <= 2 * n + 2; i += 2)
Atflo[i] = i / 2 - 1 - ell;
for (i = 2; i <= 2 * n + 2; i += 2)
Atcei[i] = i / 2 + ell;
for (i = 1; i <= 2 * n + 2; i += 2)
Atcei[i] = i / 2 + 1 + ell;
} else if (z > 0.0) {
for (i = 1; i <= 2 * n + 2; i++)
Atflo[i] = i / 2 - 1 - ell;
for (i = 2; i <= 2 * n + 2; i++)
Atcei[i] = i / 2 + ell;
Atcei[1] = 1 + ell;
} else { // z == 0
for (i = 2; i <= 2 * n + 2; i += 2)
Atflo[i] = i / 2 - 1 - ell;
for (i = 1; i <= 2 * n + 2; i += 2)
Atflo[i] = i / 2 - ell;
for (i = 2; i <= 2 * n + 2; i += 2)
Atcei[i] = i / 2 - 1 + ell;
for (i = 1; i <= 2 * n + 2; i += 2)
Atcei[i] = i / 2 + ell;
}
if (w < z)
z = w;
A[0] = A[1] = 0;
A[2] = z;
A[3] = 1 - A[2];
for (i = 4; i <= 2 * n + 1; i++)
A[i] = A[i - 2] + 1;
A[2 * n + 2] = n;
}
//========================================================================
private static double Pomeranz (int n, double x)
{
// The Pomeranz algorithm to compute the KS distribution
final double EPS = 1.0e-15;
final int ENO = 350;
final double RENO = Math.scalb (1.0, ENO); // for renormalization of V
int coreno; // counter: how many renormalizations
final double t = n * x;
double w, sum, minsum;
int i, j, k, s;
int r1, r2; // Indices i and i-1 for V[i][]
int jlow, jup, klow, kup, kup0;
double[] A = new double[2 * n + 3];
double[] Atflo = new double[2 * n + 3];
double[] Atcei = new double[2 * n + 3];
double[][] V = new double[2][n + 2];
double[][] H = new double[4][n + 2]; // = pow(w, j) / Factorial(j)
CalcFloorCeil (n, t, A, Atflo, Atcei);
for (j = 1; j <= n + 1; j++)
V[0][j] = 0;
for (j = 2; j <= n + 1; j++)
V[1][j] = 0;
V[1][1] = RENO;
coreno = 1;
// Precompute H[][] = (A[j] - A[j-1]^k / k!
H[0][0] = 1;
w = 2.0 * A[2] / n;
for (j = 1; j <= n + 1; j++)
H[0][j] = w * H[0][j - 1] / j;
H[1][0] = 1;
w = (1.0 - 2.0 * A[2]) / n;
for (j = 1; j <= n + 1; j++)
H[1][j] = w * H[1][j - 1] / j;
H[2][0] = 1;
w = A[2] / n;
for (j = 1; j <= n + 1; j++)
H[2][j] = w * H[2][j - 1] / j;
H[3][0] = 1;
for (j = 1; j <= n + 1; j++)
H[3][j] = 0;
r1 = 0;
r2 = 1;
for (i = 2; i <= 2 * n + 2; i++) {
jlow = (int) (2 + Atflo[i]);
if (jlow < 1)
jlow = 1;
jup = (int) (Atcei[i]);
if (jup > n + 1)
jup = n + 1;
klow = (int) (2 + Atflo[i - 1]);
if (klow < 1)
klow = 1;
kup0 = (int) (Atcei[i - 1]);
// Find to which case it corresponds
w = (A[i] - A[i - 1]) / n;
s = -1;
for (j = 0; j < 4; j++) {
if (Math.abs (w - H[j][1]) <= EPS) {
s = j;
break;
}
}
minsum = RENO;
r1 = (r1 + 1) & 1; // i - 1
r2 = (r2 + 1) & 1; // i
for (j = jlow; j <= jup; j++) {
kup = kup0;
if (kup > j)
kup = j;
sum = 0;
for (k = kup; k >= klow; k--)
sum += V[r1][k] * H[s][j - k];
V[r2][j] = sum;
if (sum < minsum)
minsum = sum;
}
if (minsum < 1.0e-280) {
// V is too small: renormalize to avoid underflow of probabilities
for (j = jlow; j <= jup; j++)
V[r2][j] *= RENO;
coreno++; // keep track of log of RENO
}
}
sum = V[r2][n + 1];
w = getLogFactorial (n) - coreno * ENO * num_Ln2 + Math.log (sum);
if (w >= 0.)
return 1.;
return Math.exp (w);
}
//========================================================================
private static double cdfSpecial (int n, double x)
{
// The KS distribution is known exactly for these cases
// For nx^2 > 18, fbar(n, x) is smaller than 5e-16
if ((n * x * x >= 18.0) || (x >= 1.0))
return 1.0;
if (x <= 0.5 / n)
return 0.0;
if (n == 1)
return 2.0 * x - 1.0;
if (x <= 1.0 / n) {
double t = 2.0 * x - 1.0 / n;
double w;
if (n <= NFACT) {
w = Factorial[n];
return w * Math.pow (t, (double) n);
}
w = getLogFactorial (n) + n * Math.log (t);
return Math.exp (w);
}
if (x >= 1.0 - 1.0 / n) {
return 1.0 - 2.0 * Math.pow (1.0 - x, (double) n);
}
return -1.0;
}
//========================================================================
/**
* Computes the cumulative probability P[D_n <= x] of the
* Kolmogorov-Smirnov distribution with sample size n at x.
* It returns at least 13 decimal digits of precision for n <= 140,
* at least 5 decimal digits of precision for 140 < n <= 100000,
* and a few correct decimal digits for n > 100000.
*
* @param n sample size
* @param x value of Kolmogorov-Smirnov statistic
* @return cumulative probability
*/
public static double cdf (int n, double x)
{
double u = cdfSpecial (n, x);
if (u >= 0.0)
return u;
final double w = n * x * x;
if (n <= NEXACT) {
if (w < 0.754693)
return DurbinMatrix (n, x);
if (w < 4.0)
return Pomeranz (n, x);
return 1.0 - fbar (n, x);
}
// if (n * x * sqrt(x) <= 1.4)
if ((w * x * n <= 2.0) && (n <= NKOLMO))
return DurbinMatrix(n, x);
return Pelz (n, x);
}
//=========================================================================
private static double fbarSpecial (int n, double x)
{
final double w = n * x * x;
if ((w >= 370.0) || (x >= 1.0))
return 0.0;
if ((w <= 0.0274) || (x <= 0.5 / n))
return 1.0;
if (n == 1)
return 2.0 - 2.0 * x;
if (x <= 1.0 / n) {
double v;
double t = 2.0 * x - 1.0 / n;
if (n <= NFACT) {
v = Factorial[n];
return 1.0 - v * Math.pow (t, (double) n);
}
v = getLogFactorial (n) + n * Math.log (t);
return 1.0 - Math.exp (v);
}
if (x >= 1.0 - 1.0 / n) {
return 2.0 * Math.pow (1.0 - x, (double) n);
}
return -1.0;
}
//========================================================================
/** Computes the complementary cumulative probability P[D_n >= x] of the
* Kolmogorov-Smirnov distribution with sample size n at x.
* It returns at least 10 decimal digits of precision for n <= 140,
* at least 5 decimal digits of precision for 140 < n <= 200000,
* and a few correct decimal digits for n > 200000.
*
* @param n sample size
* @param x value of Kolmogorov-Smirnov statistic
* @return complementary cumulative probability
*/
public static double fbar (int n, double x)
{
double v = fbarSpecial (n, x);
if (v >= 0.0)
return v;
final double w = n * x * x;
if (n <= NEXACT) {
if (w < 4.0)
return 1.0 - cdf (n, x);
else
return 2.0 * KSPlusbarUpper (n, x);
}
if (w >= 2.2)
return 2.0 * KSPlusbarUpper (n, x);
return 1.0 - cdf (n, x);
}
/*=========================================================================
The following implements the Durbin matrix algorithm and was programmed
by G. Marsaglia, Wai Wan Tsang and Jingbo Wong in C.
I have translated their program in Java. I have made small modifications
in their program. (Richard Simard)
=========================================================================*/
/*
The C program to compute Kolmogorov's distribution
K(n,d) = Prob(D_n < d), where
D_n = max(x_1-0/n,x_2-1/n...,x_n-(n-1)/n,1/n-x_1,2/n-x_2,...,n/n-x_n)
with x_17 digit accuracy in the right tail
if (false) {
s = d * d * n;
if (s > 7.24 || (s > 3.76 && n > 99))
return 1 - 2 * Math.exp (-(2.000071 + .331 / Math.sqrt (n) +
1.409 / n) * s);
}
k = (int) (n * d) + 1;
m = 2 * k - 1;
h = k - n * d;
H = new double[m * m];
Q = new double[m * m];
pQ = new int[1];
for (i = 0; i < m; i++)
for (j = 0; j < m; j++)
if (i - j + 1 < 0)
H[i * m + j] = 0;
else
H[i * m + j] = 1;
for (i = 0; i < m; i++) {
H[i * m] -= Math.pow (h, (double)(i + 1));
H[(m - 1) * m + i] -= Math.pow (h, (double)(m - i));
}
H[(m - 1) * m] += (2 * h - 1 > 0 ? Math.pow (2 * h - 1, (double) m) : 0);
for (i = 0; i < m; i++)
for (j = 0; j < m; j++)
if (i - j + 1 > 0)
for (g = 1; g <= i - j + 1; g++)
H[i * m + j] /= g;
eH = 0;
mPower (H, eH, Q, pQ, m, n);
s = Q[(k - 1) * m + k - 1];
for (i = 1; i <= n; i++) {
s = s * (double) i / n;
if (s < INORM) {
s *= NORM;
pQ[0] -= LOGNORM;
}
}
s *= Math.pow (10., (double) pQ[0]);
return s;
}
private static void mMultiply (double[] A, double[] B, double[] C, int m)
{
int i, j, k;
double s;
for (i = 0; i < m; i++)
for (j = 0; j < m; j++) {
s = 0.;
for (k = 0; k < m; k++)
s += A[i * m + k] * B[k * m + j];
C[i * m + j] = s;
}
}
private static void renormalize (double[] V, int m, int[] p)
{
for (int i = 0; i < m * m; i++)
V[i] *= INORM;
p[0] += LOGNORM;
}
private static void mPower (double[] A, int eA, double[] V, int[] eV, int m,
int n)
{
int i;
if (n == 1) {
for (i = 0; i < m * m; i++)
V[i] = A[i];
eV[0] = eA;
return;
}
mPower (A, eA, V, eV, m, n / 2);
double[] B = new double[m * m];
int[] pB = new int[1];
mMultiply (V, V, B, m);
pB[0] = 2 * (eV[0]);
if (B[(m / 2) * m + (m / 2)] > NORM)
renormalize (B, m, pB);
if (n % 2 == 0) {
for (i = 0; i < m * m; i++)
V[i] = B[i];
eV[0] = pB[0];
} else {
mMultiply (A, B, V, m);
eV[0] = eA + pB[0];
}
if (V[(m / 2) * m + (m / 2)] > NORM)
renormalize (V, m, eV);
}
//======================================================================
public static void main (String[] args) {
double x, y, z;
final int K = 100;
int n = 60;
System.out.printf("n = %d%n%n", n);
System.out.println(" x cdf fbar");
for (int j = 0; j <= K; j++) {
x = (double) j / K;
y = cdf (n, x);
z = fbar (n, x);
System.out.printf ("%8.3f %22.15g %22.15g%n", x, y, z);
}
}
}