/******************************************************************** * * File: KolmogorovSmirnovDist.c * Environment: ISO C99 or ANSI C89 * Author: Richard Simard * Organization: DIRO, Université de Montréal * Date: 1 February 2012 * Version 1.1 * 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 . =====================================================================*/ #include "KolmogorovSmirnovDist.h" #include #include #define num_Pi 3.14159265358979323846 /* PI */ #define num_Ln2 0.69314718055994530941 /* log(2) */ /* 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. */ #define NEXACT 500 #define NKOLMO 100000 /* The Durbin matrix algorithm for the Kolmogorov-Smirnov distribution */ static double DurbinMatrix (int n, double d); /*========================================================================*/ #if 0 /* For ANSI C89 only, not for ISO C99 */ #define MAXI 50 #define EPSILON 1.0e-15 double log1p (double x) { /* returns a value equivalent to log(1 + x) accurate also for small x. */ if (fabs (x) > 0.1) { return log (1.0 + x); } else { double term = x; double sum = x; int s = 2; while ((fabs (term) > EPSILON * fabs (sum)) && (s < MAXI)) { term *= -x; sum += term / s; s++; } return sum; } } #undef MAXI #undef EPSILON #endif /*========================================================================*/ #define MFACT 30 /* The natural logarithm of factorial n! for 0 <= n <= MFACT */ static double LnFactorial[MFACT + 1] = { 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 }; /*------------------------------------------------------------------------*/ 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) * log (x) - x) + 9.1893853320467E-1 + z / x; return z; } } /*------------------------------------------------------------------------*/ static double rapfac (int n) { /* Computes n! / n^n */ int i; double res = 1.0 / n; for (i = 2; i <= n; i++) { res *= (double) i / n; } return res; } /*========================================================================*/ static double **CreateMatrixD (int N, int M) { int i; double **T2; T2 = (double **) malloc (N * sizeof (double *)); T2[0] = (double *) malloc (N * M * sizeof (double)); for (i = 1; i < N; i++) T2[i] = T2[0] + i * M; return T2; } static void DeleteMatrixD (double **T) { free (T[0]); free (T); } /*========================================================================*/ 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 * exp (-z); if (v >= 1.0) return 1.0; return v; } /*-------------------------------------------------------------------------*/ static double KSPlusbarUpper (int n, double x) { /* Compute the probability of the KS+ distribution in the upper tail using Smirnov's stable formula */ const double EPSILON = 1.0E-12; double q; double Sum = 0.0; double term; double t; double LogCom; double LOGJMAX; int j; int jdiv; int jmax = (int) (n * (1.0 - x)); if (n > 200000) return KSPlusbarAsymp (n, x); /* Avoid log(0) for j = jmax and q ~ 1.0 */ if ((1.0 - x - (double) jmax / n) <= 0.0) jmax--; if (n > 3000) jdiv = 2; else jdiv = 3; j = jmax / jdiv + 1; LogCom = getLogFactorial (n) - getLogFactorial (j) - getLogFactorial (n - j); LOGJMAX = LogCom; while (j <= jmax) { q = (double) j / n + x; term = LogCom + (j - 1) * log (q) + (n - j) * log1p (-q); t = exp (term); Sum += t; LogCom += log ((double) (n - j) / (j + 1)); if (t <= Sum * EPSILON) break; j++; } j = jmax / jdiv; LogCom = LOGJMAX + log ((double) (j + 1) / (n - j)); while (j > 0) { q = (double) j / n + x; term = LogCom + (j - 1) * log (q) + (n - j) * log1p (-q); t = exp (term); Sum += t; LogCom += log ((double) j / (n - j + 1)); if (t <= Sum * EPSILON) break; j--; } Sum *= x; /* add the term j = 0 */ Sum += exp (n * log1p (-x)); return Sum; } /*========================================================================*/ 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 */ const int JMAX = 20; const double EPS = 1.0e-10; const double C = 2.506628274631001; /* sqrt(2*Pi) */ const double C2 = 1.2533141373155001; /* sqrt(Pi/2) */ const double PI2 = num_Pi * num_Pi; const double PI4 = PI2 * PI2; const double RACN = sqrt ((double) n); const double z = RACN * x; const double z2 = z * z; const double z4 = z2 * z2; const double z6 = z4 * z2; const 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 = exp (-ti * ti * w); sum += term; j++; } sum *= C / z; term = 1; tom = 0; j = 0; while (j <= JMAX && fabs (term) > EPS * fabs (tom)) { ti = j + 0.5; term = (PI2 * ti * ti - z2) * exp (-ti * ti * w); tom += term; j++; } sum += tom * C2 / (RACN * 3.0 * z4); term = 1; tom = 0; j = 0; while (j <= JMAX && fabs (term) > EPS * fabs (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 *= 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 * exp (-ti * ti * w); tom += term; j++; } sum -= tom * C2 / (n * 18.0 * z * z2); term = 1; tom = 0; j = 0; while (j <= JMAX && fabs (term) > EPS * fabs (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 *= exp (-ti * w); tom += term; j++; } sum += tom * C2 / (RACN * n * 3240.0 * z4 * z6); term = 1; tom = 0; j = 1; while (j <= JMAX && fabs (term) > EPS * fabs (tom)) { ti = j * j; term = (3 * PI2 * ti * z2 - PI4 * ti * ti) * exp (-ti * w); tom += term; j++; } sum += tom * C2 / (RACN * n * 108.0 * z6); return sum; } /*=========================================================================*/ 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 = 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; } /*========================================================================*/ static double Pomeranz (int n, double x) { /* The Pomeranz algorithm to compute the KS distribution */ const double EPS = 1.0e-15; const int ENO = 350; const double RENO = ldexp (1.0, ENO); /* for renormalization of V */ int coreno; /* counter: how many renormalizations */ const 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; double *Atflo; double *Atcei; double **V; double **H; /* = pow(w, j) / Factorial(j) */ A = (double *) calloc ((size_t) (2 * n + 3), sizeof (double)); Atflo = (double *) calloc ((size_t) (2 * n + 3), sizeof (double)); Atcei = (double *) calloc ((size_t) (2 * n + 3), sizeof (double)); V = (double **) CreateMatrixD (2, n + 2); H = (double **) CreateMatrixD (4, n + 2); 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! for speed */ 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 = 2 + (int) Atflo[i]; if (jlow < 1) jlow = 1; jup = (int) Atcei[i]; if (jup > n + 1) jup = n + 1; klow = 2 + (int) 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 (fabs (w - H[j][1]) <= EPS) { s = j; break; } } /* assert (s >= 0, "Pomeranz: s < 0"); */ 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]; free (A); free (Atflo); free (Atcei); DeleteMatrixD (H); DeleteMatrixD (V); w = getLogFactorial (n) - coreno * ENO * num_Ln2 + log (sum); if (w >= 0.) return 1.; return exp (w); } /*========================================================================*/ static double cdfSpecial (int n, double x) { /* The KS distribution is known exactly for these cases */ /* For nx^2 > 18, KSfbar(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 * n - 1.0; double w; if (n <= NEXACT) { w = rapfac (n); return w * pow (t, (double) n); } w = getLogFactorial (n) + n * log (t / n); return exp (w); } if (x >= 1.0 - 1.0 / n) { return 1.0 - 2.0 * pow (1.0 - x, (double) n); } return -1.0; } /*========================================================================*/ double KScdf (int n, double x) { const double w = n * x * x; double u = cdfSpecial (n, x); if (u >= 0.0) return u; if (n <= NEXACT) { if (w < 0.754693) return DurbinMatrix (n, x); if (w < 4.0) return Pomeranz (n, x); return 1.0 - KSfbar (n, x); } if ((w * x * n <= 7.0) && (n <= NKOLMO)) return DurbinMatrix (n, x); return Pelz (n, x); } /*=========================================================================*/ static double fbarSpecial (int n, double x) { const 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 z; double t = 2.0 * x * n - 1.0; if (n <= NEXACT) { z = rapfac (n); return 1.0 - z * pow (t, (double) n); } z = getLogFactorial (n) + n * log (t / n); return 1.0 - exp (z); } if (x >= 1.0 - 1.0 / n) { return 2.0 * pow (1.0 - x, (double) n); } return -1.0; } /*========================================================================*/ double KSfbar (int n, double x) { const double w = n * x * x; double v = fbarSpecial (n, x); if (v >= 0.0) return v; if (n <= NEXACT) { if (w < 4.0) return 1.0 - KScdf (n, x); else return 2.0 * KSPlusbarUpper (n, x); } if (w >= 2.65) return 2.0 * KSPlusbarUpper (n, x); return 1.0 - KScdf (n, x); } /*========================================================================= The following implements the Durbin matrix algorithm and was programmed by G. Marsaglia, Wai Wan Tsang and Jingbo Wong. 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 0 s = d * d * n; if (s > 7.24 || (s > 3.76 && n > 99)) return 1 - 2 * exp (-(2.000071 + .331 / sqrt (n) + 1.409 / n) * s); #endif k = (int) (n * d) + 1; m = 2 * k - 1; h = k - n * d; H = (double *) malloc ((m * m) * sizeof (double)); Q = (double *) malloc ((m * m) * sizeof (double)); 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] -= pow (h, (double) (i + 1)); H[(m - 1) * m + i] -= pow (h, (double) (m - i)); } H[(m - 1) * m] += (2 * h - 1 > 0 ? 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, &eQ, 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; eQ -= LOGNORM; } } s *= pow (10., (double) eQ); free (H); free (Q); return s; } 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; } } static void renormalize (double *V, int m, int *p) { int i; for (i = 0; i < m * m; i++) V[i] *= INORM; *p += LOGNORM; } static void mPower (double *A, int eA, double *V, int *eV, int m, int n) { double *B; int eB, i; if (n == 1) { for (i = 0; i < m * m; i++) V[i] = A[i]; *eV = eA; return; } mPower (A, eA, V, eV, m, n / 2); B = (double *) malloc ((m * m) * sizeof (double)); mMultiply (V, V, B, m); eB = 2 * (*eV); if (B[(m / 2) * m + (m / 2)] > NORM) renormalize (B, m, &eB); if (n % 2 == 0) { for (i = 0; i < m * m; i++) V[i] = B[i]; *eV = eB; } else { mMultiply (A, B, V, m); *eV = eA + eB; } if (V[(m / 2) * m + (m / 2)] > NORM) renormalize (V, m, eV); free (B); } /*=========================================================================*/ #if 0 #include int main (void) { double x, y, z; const int K = 100; int n = 60; int j; printf ("n = %5d\n\n", n); printf (" x cdf fbar\n"); for (j = 0; j <= K; j++) { x = (double) j / K; y = KScdf (n, x); z = KSfbar (n, x); printf ("%8.3g %22.15g %22.15g\n", x, y, z); } return 0; } #endif