/* *
* 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 ContinuousDistribution} for the
* Kolmogorov-Smirnov distribution with parameter n.
* Given an empirical distribution Fn with n independent observations and
* a continuous distribution F(x), the two-sided statistic is defined as
*
*
*
* Dn = sup-∞ <= x <= ∞| Fn(x) - F(x)| = {Dn+, Dn-},
*
* where Dn+ and Dn- are the + and - statistics as defined in
* equations and on page ![[*]](file:/usr/share/latex2html/icons/crossref.png)
of this guide.
* This class implements a high precision version of the
* distribution
* P[Dn <= x]; it is a Java translation of the C program
* written in. According to its authors, it should give
* 13 decimal digits of precision. It is extremely slow
* for large values of n.
*
*/
public class KolmogorovSmirnovDist extends ContinuousDistribution {
protected int n;
protected static final int NLIM = 20;
// 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;
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.
* Restriction: n >= 1.
*
*/
public KolmogorovSmirnovDist (int n) {
setN (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);
}
private static double dclem (int n, double x, double EPS) {
return (cdf(n, x + EPS) - cdf(n, x - EPS)) / (2.0 * EPS);
}
protected static double densConnue (int n, double x) {
if ((x >= 1.0) || (x <= 0.5 / n))
return 0.0;
if (n == 1)
return 2.0;
if (x <= 1.0 / n) {
double w;
double t = 2.0 * x - 1.0 / n;
if (n <= NLIM) {
w = 2.0 * n * Num.factorial (n);
w *= Math.pow (t, (double) (n - 1));
return w;
}
w = Num.lnFactorial (n) + (n-1) * Math.log (t);
return 2*n*Math.exp (w);
}
if (x >= 1.0 - 1.0 / n)
return 2.0 * n * Math.pow (1.0 - x, (double) (n - 1));
return -1.0;
}
/**
* 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;
double EPS = 1.0 / 200.0;
final double D1 = dclem(n, x, EPS);
final double D2 = dclem(n, x, 2.0 * EPS);
Res = D1 + (D1 - D2) / 3.0;
if (Res <= 0.0)
return 0.0;
return Res;
}
/*=========================================================================
The following implements the Durbin matrix algorithm and was programmed
in C by G. Marsaglia, Wai Wan Tsang and Jingbo Wong in C.
I have translated their program in Java. Only small modifications have
been made in their program; the most important is to prevent the return
of NAN or infinite values in some regions. (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_1 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);
}
protected static double DurbinMatrix (int n, double d)
{
int k, m, i, j, g, eH;
double h, s;
double[] H;
double[] Q;
int[] pQ;
//Omit next two lines if you require >7 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;
}
protected static double cdfConnu (int n, double x) {
// For nx^2 > 18, barF(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 w;
double t = 2.0 * x - 1.0 / n;
if (n <= NLIM) {
w = Num.factorial (n);
return w * Math.pow (t, (double) n);
}
w = Num.lnFactorial (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 distribution function F(x) with parameter n
* using Durbin's matrix formula. It is a translation of the
* C program in;
* according to its authors, it returns 13 decimal digits
* of precision. It is extremely slow for large n.
*
*/
public static double cdf (int n, double x) {
double Res = cdfConnu(n,x);
if (Res != -1.0)
return Res;
return DurbinMatrix (n, x);
}
protected static double barFConnu (int n, double x) {
if ((n * x * x >= 370.0) || (x >= 1.0))
return 0.0;
if ((n * x * x <= 0.0274) || (x <= 0.5 / n))
return 1.0;
if (n == 1)
return 2.0 - 2.0 * x;
if (x <= 1.0 / n) {
double w;
final double t = 2.0 * x - 1.0 / n;
if (n <= NLIM) {
w = Num.factorial (n);
return 1.0 - w * Math.pow (t, (double) n);
}
w = Num.lnFactorial (n) + n * Math.log (t);
return 1.0 - Math.exp (w);
}
if (x >= 1.0 - 1.0 / n) {
return 2.0 * Math.pow (1.0 - x, (double) n);
}
return -1.0;
}
/**
* Computes the complementary distribution function bar(F)(x)
* with parameter n. Simply returns 1 - cdf(n,x). It is not precise in
* the upper tail.
*
*/
public static double barF (int n, double x) {
double h = barFConnu(n, x);
if (h >= 0.0)
return h;
h = 1.0 - cdf(n, x);
if (h >= 0.0)
return h;
return 0.0;
}
protected static double inverseConnue (int n, double u) {
if (n <= 0)
throw new IllegalArgumentException ("n <= 0");
if (u < 0.0 || u > 1.0)
throw new IllegalArgumentException ("u must be in [0,1]");
if (u == 1.0)
return 1.0;
if (u == 0.0)
return 0.5/n;
if (n == 1)
return (u + 1.0) / 2.0;
final double NLNN = n*Math.log (n);
final double LNU = Math.log(u) - Num.lnFactorial (n);
if (LNU <= -NLNN){
double t = 1.0/n*(LNU);
return 0.5 * (Math.exp(t) + 1.0/n);
}
if (u >= 1.0 - 2.0 / Math.exp (NLNN))
return 1.0 - Math.pow((1.0-u)/2.0, 1.0/n);
return -1.0;
}
/**
* 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-10);
}
/**
* Returns the parameter n of this object.
*
*/
public int getN() {
return n;
}
/**
* Sets the parameter n of this object.
*
*/
public void setN (int n) {
if (n <= 0)
throw new IllegalArgumentException ("n <= 0");
this.n = n;
supportA = 0.5 / n;
supportB = 1.0;
}
/**
* Returns an array containing the parameter n of this object.
*
*
*/
public double[] getParams () {
double[] retour = {n};
return retour;
}
public String toString () {
return getClass().getSimpleName() + " : n = " + n;
}
}