/* *
* Title:
* Description:
* Copyright: Copyright (c)
* Company:
* @author
* @version 1.0
*/
package umontreal.iro.lecuyer.probdist;
import umontreal.iro.lecuyer.util.*;
import umontreal.iro.lecuyer.functions.MathFunction;
/**
* Extends the class {@link KolmogorovSmirnovDist} for the distribution.
* The methods of this class are much faster than those of class
* {@link KolmogorovSmirnovDist}.
*
*/
public class KolmogorovSmirnovDistQuick extends KolmogorovSmirnovDist {
private static final double PI2 = Math.PI * Math.PI;
private static final double PI4 = PI2 * PI2;
private static final double C2PI = 2.506628274631001; // sqrt(2*Pi)
private static final double DPI2 = 1.2533141373155001; // sqrt(Pi/2)
/* 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;
private static class Function implements MathFunction {
protected int n;
protected double u;
public Function (int n, double u) {
this.n = n;
this.u = u;
}
public double evaluate (double x) {
return u - cdf(n,x);
}
}
/**
* Constructs a distribution with parameter n.
*
*/
public KolmogorovSmirnovDistQuick (int n) {
super (n);
}
public double density (double x) {
return density (n, x);
}
public double cdf (double x) {
return cdf (n, x);
}
public double barF (double x) {
return barF (n, x);
}
public double inverseF (double u) {
return inverseF (n, u);
}
/**
* Computes the density for the distribution with parameter n.
*
*/
public static double density (int n, double x) {
double Res = densConnue(n,x);
if (Res != -1.0)
return Res;
final double EPS = 1.0 / Num.TWOEXP[6];
Res = (cdf(n, x + EPS) - cdf(n, x - EPS)) / (2.0 * EPS);
if (Res <= 0.0)
return 0.0;
return Res;
}
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 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 *= C2PI / 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 * DPI2 / (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 * DPI2 / (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 * DPI2 / (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 * DPI2 / (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 * DPI2 / (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 = Num.lnFactorial (n) - coreno * ENO * Num.LN2 + Math.log (sum);
if (w >= 0.)
return 1.;
return Math.exp (w);
}
//========================================================================
/**
* Computes the distribution function
* u = P[Dn <= x] with
* parameter n, using the program described in.
* This method uses Pomeranz's recursion algorithm and the Durbin matrix algorithm for
* n <= 140, which returns at least 13 decimal digits of precision. It uses
* the Pelz-Good asymptotic expansion in the central part of
* the range for n > 140 and returns at least 5 decimal digits of precision
* everywhere for
* 140 < n <= 100000. For
* n > 100000, it returns
* at least 5 decimal digits of precision for all
* u > 10-16, and a
* few correct decimals when
* u <= 10-16.
* For a given n > 140, the precision increases as x increases.
* This method is much faster than method cdf of
* {@link KolmogorovSmirnovDist} for moderate or large n.
* Restriction: n >= 1.
*
*/
public static double cdf (int n, double x) {
double u = cdfConnu (n, x);
if (u >= 0.0)
return u;
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 - barF (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);
}
/**
* Computes the complementary distribution
*
* P[Dn >= x] with parameter n,
* in a form that is more precise in the upper tail,
* using the program described in.
* It returns at least 10 decimal digits of precision everywhere for all
* n <= 140,
* at least 5 decimal digits of precision for
* 140 < n <= 200000,
* and a few correct digits (1 to 5) for
* n > 200000.
* This method is much faster and more precise for x close to 1, than
* method barF of
* {@link KolmogorovSmirnovDist} for moderate or large n.
* Restriction: n >= 1.
*
*/
public static double barF (int n, double x) {
double v = barFConnu (n, x);
if (v >= 0.0)
return v;
double w = n * x * x;
if (n <= NEXACT) {
if (w < 4.0)
return 1.0 - cdf (n, x);
else
return 2.0 * KolmogorovSmirnovPlusDist.KSPlusbarUpper(n, x);
}
if (w >= 2.2) {
if (n <= 200000)
return 2.0 * KolmogorovSmirnovPlusDist.KSPlusbarUpper (n, x);
return 2.0 * KolmogorovSmirnovPlusDist.KSPlusbarAsymp (n, x);
}
return 1.0 - cdf (n, x);
}
/**
* Computes the inverse
* x = F-1(u) of the
* distribution F(x) with parameter n.
*
*/
public static double inverseF (int n, double u) {
double Res = inverseConnue(n,u);
if (Res != -1.0)
return Res;
Function f = new Function (n,u);
return RootFinder.brentDekker (0.5/n, 1.0, f, 1e-5);
}
}
Ads help maintain this website.