/* *
* 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
* Watson U distribution (see).
* Given a sample of n independent uniforms ui over [0, 1],
* the Watson statistic Un2 is defined by
*
*
*
* Wn2 = 1/12n + ∑j=1n{u(j) - (j-0.5)/n}2,
*
Un2 = Wn2 - n(bar(u)n -1/2)2.
*
* where the u(j) are the ui sorted in increasing order, and
* bar(u)n is the average of the observations ui.
* The distribution function (the cumulative probabilities)
* is defined as
* Fn(x) = P[Un2 <= x].
*
*/
public class WatsonUDist extends ContinuousDistribution {
private static final double XSEPARE = 0.15;
private static final double PI = Math.PI;
private static final int JMAX = 10;
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);
}
}
private static double cdfn (int n, double x) {
// The 1/n correction for the cdf, for x < XSEPARE
double terme;
double v = Math.exp (-0.125/x);
double somme = 0;
int j = 0;
do {
double a = (2*j + 1)*(2*j + 1);
terme = Math.pow (v, (double)(2*j + 1)*(2*j + 1));
double der = terme*(a - 4.0*x)/(8.0*x*x);
somme += (5.0*x - 1.0/12.0) * der / 12.0;
der = terme* (a*a - 24.0*a*x + 48.0*x*x)/ (64.0*x*x*x*x);
somme += x*x*der/6.0;
++j;
} while (!(terme <= Math.abs(somme) * Num.DBL_EPSILON || j > JMAX));
if (j > JMAX)
System.err.println (x + ": watsonU: somme 1/n has not converged");
v = -2.0*somme/(n*Math.sqrt (2.0*PI*x));
return v;
}
/**
* Constructs a Watson U distribution for a sample of size n.
*
*/
public WatsonUDist (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 WatsonUDist.getMean (n);
}
public double getVariance() {
return WatsonUDist.getVariance (n);
}
public double getStandardDeviation() {
return WatsonUDist.getStandardDeviation (n);
}
/**
* Computes the density of the Watson U distribution
* with parameter n.
*
*/
public static double density (int n, double x) {
if (n < 2)
throw new IllegalArgumentException ("n < 2");
if (x <= 1.0/(12.0*n) || x >= n/12.0 || x >= XBIG)
return 0.0;
final double EPS = 1.0 / 100.0;
return (cdf(n, x + EPS) - cdf(n, x - EPS)) / (2.0 * EPS);
/*
// This is the asymptotic density for n -> infinity
int j;
double signe;
double v;
double terme;
double somme;
if (x > XSEPARE) {
// this series converges rapidly for x > 0.15
v = Math.exp (-(x*2.0*PI*PI));
signe = 1.0;
somme = 0.0;
j = 1;
do {
terme = j*j*Math.pow (v, (double)j*j);
somme += signe*terme;
signe = -signe;
++j;
} while (!(terme < Num.DBL_EPSILON || j > JMAX));
if (j > JMAX)
System.err.println ("watsonU: sum1 has not converged");
return 4.0*PI*PI*somme;
}
// this series converges rapidly for x <= 0.15
v = Math.exp (-0.125/x);
somme = v;
double somme2 = v;
j = 2;
do {
terme = Math.pow (v, (double)(2*j - 1)*(2*j - 1));
somme += terme;
somme2 += terme * (2*j - 1)*(2*j - 1);
++j;
} while (!(terme <= somme2 * Num.DBL_EPSILON || j > JMAX));
if (j > JMAX)
System.err.println ("watsonU: sum2 has not converged");
final double RACINE = Math.sqrt (2.0*PI*x);
return -somme/(x*RACINE) + 2*somme2 * 0.125/ ((x*x) * RACINE);
*/
}
/**
* Computes the Watson U distribution function, i.e. returns
*
* P[Un2 <= x], where Un2 is the Watson statistic
* defined in. We use the asymptotic distribution for
*
* n -> ∞, plus a correction in O(1/n), as given in.
*
*/
public static double cdf (int n, double x) {
if (n < 2)
throw new IllegalArgumentException ("n < 2");
if (x <= 1.0/(12.0*n))
return 0.0;
if (x > 3.95 || x >= n/12.0)
return 1.0;
if (2 == n) {
if (x <= 1.0/24.0)
return 0.0;
if (x >= 1.0/6.0)
return 1.0;
return 2.0*Math.sqrt(2.0*x - 1.0/12.0);
}
if (x > XSEPARE)
return 1.0 - barF (n, x);
// this series converges rapidly for x <= 0.15
double terme;
double v = Math.exp (-0.125/x);
double somme = v;
int j = 2;
do {
terme = Math.pow (v, (double)(2*j - 1)*(2*j - 1));
somme += terme;
++j;
} while (!(terme <= somme * Num.DBL_EPSILON || j > JMAX));
if (j > JMAX)
System.err.println (x + ": watsonU: sum2 has not converged");
v = 2.0*somme/Math.sqrt (2.0*PI*x);
v += cdfn(n, x);
if (v >= 1.0)
return 1.0;
if (v <= 0.0)
return 0.0;
return v;
}
/**
* Computes the complementary distribution function
* bar(F)n(x),
* where Fn is the Watson U distribution with parameter n.
*
*/
public static double barF (int n, double x) {
if (n < 2)
throw new IllegalArgumentException ("n < 2");
if (x <= 1.0/(12.0*n))
return 1.0;
if (x >= XBIG || x >= n/12.0)
return 0.0;
if (2 == n)
return 1.0 - 2.0*Math.sqrt(2.0*x - 1.0/12.0);
if (x > XSEPARE) {
// this series converges rapidly for x > 0.15
double terme, ter;
double v = Math.exp (-2.0*PI*PI*x);
double signe = 1.0;
double somme = 0.0;
double son = 0.0;
int j = 1;
do {
terme = Math.pow (v, (double)j*j);
somme += signe*terme;
double h = 2*j*PI*x;
ter = (5.0*x - h*h - 1.0/12.0)*j*j;
son += signe*terme*ter;
signe = -signe;
++j;
} while (!(terme < Num.DBL_EPSILON || j > JMAX));
if (j > JMAX)
System.err.println (x + ": watsonU: sum1 has not converged");
v = 2.0*somme + PI*PI*son/(3.0*n);
if (v <= 0.0)
return 0.0;
if (v >= 1.0)
return 1.0;
return v;
}
return (1.0 - cdf(n, x));
}
/**
* Computes
* x = Fn-1(u), where Fn is the
* Watson U distribution with parameter n.
*
*/
public static double inverseF (int n, double u) {
if (n < 2)
throw new IllegalArgumentException ("n < 2");
if (u < 0.0 || u > 1.0)
throw new IllegalArgumentException ("u must be in [0,1]");
if (u >= 1.0)
return n/12.0;
if (u <= 0.0)
return 1.0/(12.0*n);
if (2 == n)
return 1.0/24.0 + u*u/8.0;
Function f = new Function (n,u);
return RootFinder.brentDekker (0.0, 2.0, f, 1e-7);
}
/**
* Returns the mean of the Watson U distribution with
* parameter n.
*
* @return Returns the mean
*
*/
public static double getMean (int n) {
return 1.0/12.0;
}
/**
* Returns the variance of the Watson U distribution with
* parameter n.
*
* @return the variance
*
*/
public static double getVariance (int n) {
return (n - 1)/(360.0*n);
}
/**
* Returns the standard deviation of the Watson U
* distribution with parameter n.
*
* @return the standard deviation
*
*/
public static double getStandardDeviation (int n) {
return Math.sqrt (WatsonUDist.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 < 2)
throw new IllegalArgumentException ("n < 2");
this.n = n;
supportA = 1.0/(12.0*n);
supportB = n/12.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;
}
}
Ads help maintain this website.