umontreal.iro.lecuyer.probdist.CramerVonMisesDist 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: CramerVonMisesDist
* Description: Cramér-von Mises 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
* Cramér-von Mises distribution (see).
* Given a sample of n independent uniforms Ui over [0, 1],
* the Cramér-von Mises statistic Wn2 is defined by
*
*
*
* Wn2 = 1/12n + ∑j=1n(U(j) - (j-0.5)/n)2,
*
* where the U(j) are the Ui sorted in increasing order.
* The distribution function (the cumulative probabilities)
* is defined as
* Fn(x) = P[Wn2 <= x].
*
*/
public class CramerVonMisesDist 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 a Cramér-von Mises distribution for a sample of size n.
*
*/
public CramerVonMisesDist (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);
}
public double getMean() {
return CramerVonMisesDist.getMean (n);
}
public double getVariance() {
return CramerVonMisesDist.getVariance (n);
}
public double getStandardDeviation() {
return CramerVonMisesDist.getStandardDeviation (n);
}
/**
* Computes the density function
* for a Cramér-von Mises distribution with parameter n.
*
*/
public static double density (int n, double x) {
if (n <= 0)
throw new IllegalArgumentException ("n <= 0");
if (x <= 1.0/(12.0*n) || x >= n/3.0)
return 0.0;
if (n == 1)
return 1.0 / Math.sqrt (x - 1.0/12.0);
if (x <= 0.002 || x > 3.95)
return 0.0;
throw new UnsupportedOperationException("density not implemented.");
}
/**
* Computes the Cramér-von Mises distribution function with parameter n.
* Returns an approximation of
* P[Wn2 <= x], where Wn2 is the
* Cramér von Mises statistic (see).
* The approximation is based on the distribution function of
* W2 = limn -> ∞Wn2, which has the following series expansion derived
* by Anderson and Darling:
*
*
*
* P(W2 <= x) = 
∑j=0∞(- 1)j
(4j+1)1/2 exp{ - 
}K1/4([tex2html_wrap_indisplay244]),
*
* where
* Kν is the modified Bessel function of the
* second kind.
* To correct for the deviation between
* P(Wn2 <= x) and
* P(W2 <= x),
* we add a correction in 1/n, obtained empirically by simulation.
* For n = 10, 20, 40, the error is less than
* 0.002, 0.001, and 0.0005, respectively, while for
* n >= 100 it is less than 0.0005.
* For
* n -> ∞, we estimate that the method returns
* at least 6 decimal digits of precision.
* For n = 1, the method uses the exact distribution:
*
* P(W12 <= x) = 2(x - 1/12)1/2 for
* 1/12 <= x <= 1/3.
*
*/
public static double cdf (int n, double x) {
if (n <= 0)
throw new IllegalArgumentException ("n <= 0");
if (n == 1) {
if (x <= 1.0/12.0)
return 0.0;
if (x >= 1.0/3.0)
return 1.0;
return 2.0*Math.sqrt (x - 1.0/12.0);
}
if (x <= 1.0/(12.0*n))
return 0.0;
if (x <= (n + 3.0)/(12.0*n*n)) {
double t = Num.lnFactorial(n) - Num.lnGamma (1.0 + 0.5*n) +
0.5*n*Math.log (Math.PI*(x - 1.0/(12.0*n)));
return Math.exp(t);
}
if (x <= 0.002)
return 0.0;
if (x > 3.95 || x >= n/3.0)
return 1.0;
final double EPSILON = Num.DBL_EPSILON;
final int JMAX = 20;
int j = 0;
double Cor, Res, arg;
double termX, termS, termJ;
termX = 0.0625/x; // 1 / (16x)
Res = 0.0;
final double A[] = {
1.0,
1.11803398875,
1.125,
1.12673477358,
1.1274116945,
1.12774323743,
1.1279296875,
1.12804477649,
1.12812074678,
1.12817350091
};
do {
termJ = 4*j + 1;
arg = termJ*termJ*termX;
termS = A[j]*Math.exp (-arg)*Num.besselK025 (arg);
Res += termS;
++j;
} while (!(termS < EPSILON || j > JMAX));
if (j > JMAX)
System.err.println ("cramerVonMises: iterations have not converged");
Res /= Math.PI*Math.sqrt (x);
// Empirical correction in 1/n
if (x < 0.0092)
Cor = 0.0;
else if (x < 0.03)
Cor = -0.0121763 + x*(2.56672 - 132.571*x);
else if (x < 0.06)
Cor = 0.108688 + x*(-7.14677 + 58.0662*x);
else if (x < 0.19)
Cor = -0.0539444 + x*(-2.22024 + x*(25.0407 - 64.9233*x));
else if (x < 0.5)
Cor = -0.251455 + x*(2.46087 + x*(-8.92836 + x*(14.0988 -
x*(5.5204 + 4.61784*x))));
else if (x <= 1.1)
Cor = 0.0782122 + x*(-0.519924 + x*(1.75148 +
x*(-2.72035 + x*(1.94487 - 0.524911*x))));
else
Cor = Math.exp (-0.244889 - 4.26506*x);
Res += Cor/n;
// This empirical correction is not very precise, so ...
if (Res <= 1.0)
return Res;
else
return 1.0;
}
/**
* Computes the complementary distribution function
* bar(F)n(x)
* with parameter n.
*
*/
public static double barF (int n, double x) {
return 1.0 - cdf(n,x);
}
/**
* Computes
* x = Fn-1(u), where Fn is the
* Cramér-von Mises 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 n/3.0;
if (u <= 0.0)
return 1.0/(12.0*n);
if (n == 1)
return 1.0/12.0 + 0.25*u*u;
Function f = new Function (n,u);
return RootFinder.brentDekker (0.0, 10.0, f, 1e-6);
}
/**
* Returns the mean of the distribution with parameter n.
*
* @return the mean
*
*/
public static double getMean (int n) {
return 1.0 / 6.0;
}
/**
* Returns the variance of the distribution with parameter n.
*
* @return variance
*
*/
public static double getVariance (int n) {
return (4.0 * n - 3.0) / (180.0 * n );
}
/**
* Returns the standard deviation of the distribution with
* parameter n.
*
* @return the standard deviation
*
*/
public static double getStandardDeviation (int n) {
return Math.sqrt(getVariance(n));
}
/**
* 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 = 1.0/(12.0*n);
supportB = n/3.0;
}
/**
* Return 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;
}
}