umontreal.iro.lecuyer.probdist.KolmogorovSmirnovPlusDist Maven / Gradle / Ivy
Go to download
Show more of this group Show more artifacts with this name
Show all versions of ssj Show documentation
Show all versions of ssj Show documentation
SSJ is a Java library for stochastic simulation, developed under the direction of Pierre L'Ecuyer,
in the Département d'Informatique et de Recherche Opérationnelle (DIRO), at the Université de Montréal.
It provides facilities for generating uniform and nonuniform random variates, computing different
measures related to probability distributions, performing goodness-of-fit tests, applying quasi-Monte
Carlo methods, collecting (elementary) statistics, and programming discrete-event simulations with both
events and processes.
The newest version!
/*
* Class: KolmogorovSmirnovPlusDist
* Description: Kolmogorov-Smirnov+ distribution
* Environment: Java
* Software: SSJ
* Copyright (C) 2001 Pierre L'Ecuyer and Université de Montréal
* Organization: DIRO, Université de Montréal
* @author
* @since
* SSJ is free software: you can redistribute it and/or modify it under
* the terms of the GNU General Public License (GPL) as published by the
* Free Software Foundation, either version 3 of the License, or
* any later version.
* SSJ 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.
* A copy of the GNU General Public License is available at
GPL licence site.
*/
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 (see).
* Given a sample of n independent uniforms Ui over [0, 1],
* the Kolmogorov-Smirnov+ statistic Dn+ and the Kolmogorov-Smirnov- statistic Dn-,
* are defined by
*
*
*
*
* Dn+
* =
* max1 <= j <= n(j/n - U(j)),
*
* Dn-
* =
* max1 <= j <= n(U(j) - (j - 1)/n),
*
*
*
*
* where the U(j) are the Ui sorted in increasing order. Both statistics
* follows the same distribution function, i.e.
*
* Fn(x) = P[Dn+ <= x] = P[Dn- <= x].
*
*/
public class KolmogorovSmirnovPlusDist extends ContinuousDistribution {
protected int n;
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 an Kolmogorov-Smirnov+ distribution for a sample of size n.
*
*/
public KolmogorovSmirnovPlusDist (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);
}
/**
* Computes the density of the Kolmogorov-Smirnov+ distribution with parameter n.
*
*/
public static double density (int n, double x) {
if (n <= 0)
throw new IllegalArgumentException (
"Calling kolmogorovSmirnovPlus with n < 1");
if (x <= 0.0 || x >= 1.0)
return 0.0;
if (n == 1)
return 1.0;
final double EPS = 1.0 / 100.0;
final double D1 = dclem(n, x, EPS);
final double D2 = dclem(n, x, 2.0 * EPS);
final double RES = D1 + (D1 - D2) / 3.0;
if (RES < 0.0)
return 0.0;
return RES;
}
/**
* Computes the Kolmogorov-Smirnov+ distribution function Fn(x) with parameter n.
* The relative error on
* Fn(x) = P[Dn+ <= x] is always less than
* 10-5.
*/
public static double cdf (int n, double x) {
if (n <= 0)
throw new IllegalArgumentException (
"Calling kolmogorovSmirnovPlus with n < 1");
if (x <= 0.0)
return 0.0;
if ((x >= 1.0) || (n*x*x >= 25.0))
return 1.0;
if (n == 1)
return x;
final double NXPARAM = 6.5; // frontier: alternating series
final int NPARAM = 4000; // frontier: non-alternating series
double q;
double Sum = 0.0;
double term;
double Njreal;
double jreal;
double LogCom = Math.log ((double)n);
int j;
int jmax;
//--------------------------------------------------------------
// the alternating series is stable and fast for n*x very small
//--------------------------------------------------------------
if (n*x <= NXPARAM) {
final double EPSILON = Num.DBL_MIN;
jmax = (int)(n*x);
int Sign = -1;
for (j = 1; j <= jmax; j++) {
jreal = j;
Njreal = n - j;
q = jreal/n - x;
// we must avoid log (0.0) for j = jmax and n*x near an integer
if (-q > EPSILON) {
term = LogCom + jreal*Math.log(-q) + (Njreal - 1.0)
*Math.log1p (-q);
Sum += Sign*Math.exp (term);
}
Sign = -Sign;
LogCom += Math.log (Njreal/(j + 1));
}
// add the term j = 0
Sum += Math.exp ((n - 1)*Math.log1p (x));
return Sum*x;
}
//-----------------------------------------------------------
// For nx > NxParam, we use the other exact series for small
// n, and the asymptotic form for n larger than NPARAM
//-----------------------------------------------------------
if (n <= NPARAM) {
jmax = (int)(n*(1.0 - x));
if (1.0 - x - (double) jmax/n <= 0.0)
--jmax;
for (j = 1; j <= jmax; j++) {
jreal = j;
Njreal = n - j;
q = jreal/n + x;
term = LogCom+(jreal - 1.0)*Math.log(q) + Njreal*Math.log1p (-q);
Sum += Math.exp (term);
LogCom += Math.log (Njreal/(jreal + 1.0));
}
Sum *= x;
// add the term j = 0; avoid log (0.0)
if (1.0 > x)
Sum += Math.exp (n*Math.log1p (-x));
return 1.0 - Sum;
}
//--------------------------
// Use an asymptotic formula
//--------------------------
term = 2.0/3.0;
q = x*x*n;
Sum = 1.0 - Math.exp (-2.0*q)*(1.0 - term*x*(1.0 - x*(1.0 - term*q)
- term/n*(0.2 - 19.0/15.0*q + term*q*q)));
return Sum;
}
private static double KSPlusbarAsymp (int n, double x) {
/* Compute the probability of the complementary KSPlus 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;
}
//-------------------------------------------------------------------------
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 = Num.lnFactorial(n) - Num.lnFactorial(j)
- Num.lnFactorial(n-j);
final double LOGJM = 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 = LOGJM + 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;
}
/**
* Computes the complementary distribution function
* bar(F)n(x)
* with parameter n.
*
*/
public static double barF (int n, double x) {
if (n <= 0)
throw new IllegalArgumentException (
"Calling kolmogorovSmirnovPlus with n < 1");
if (x <= 0.0)
return 1.0;
if ((x >= 1.0) || (n*x*x >= 365.0))
return 0.0;
if (n == 1)
return 1.0 - x;
final double NXPARAM = 6.5; // frontier: alternating series
final int NPARAM = 4000; // frontier: non-alternating series
final int NASYMP = 200000; // frontier: asymptotic
// the alternating series is stable and fast for n*x very small
if (n*x <= NXPARAM)
return 1.0 - cdf(n, x);
if (n >= NASYMP)
return KSPlusbarAsymp (n, x);
if ((n <= NPARAM) || (n*x*x > 1.0))
return KSPlusbarUpper(n,x);
return KSPlusbarAsymp (n, x);
// return (1.0 - 2.0*x/3.0)*Math.exp(-2.0*n*x*x);
}
/**
* Computes the inverse
* x = F-1(u) of the distribution with parameter n.
*
*/
public static double inverseF (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.0;
Function f = new Function (n,u);
return RootFinder.brentDekker (0.0, 1.0, f, 1e-8);
}
/**
* 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.0;
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;
}
}
© 2015 - 2025 Weber Informatics LLC | Privacy Policy