package umontreal.iro.lecuyer.gof;
import umontreal.iro.lecuyer.probdist.*;
/**
* This class is similar to {@link FDist}, except that it provides static methods
* to compute or approximate the complementary distribution function of X,
* which we define as
* bar(F)(x) = P[X >= x], instead of
* F(x) = P[X <= x].
* Note that with our definition of bar(F), one has
*
* bar(F)(x) = 1 - F(x) for continuous distributions and
*
* bar(F)(x) = 1 - F(x - 1) for discrete distributions over the integers.
*
*/
public class FBar {
private FBar() {}
private static final double EPSILONSCAN = 1.0E-7;
private static double scanGlaz (int n, double d, int m) {
int j, jmoy;
double temp;
double jr, jm1r, nr = n;
int signe;
double q = 1.0 - d;
double q4, q3, q2, q1;
double bin, binMoy;
jmoy = (int)((n + 1)*d); // max term of the Binomial
if (jmoy < m - 1)
jmoy = m - 1;
/*---------------------------------------------------------
* Compute q1: formula (2.5) in Glaz (1989)
* Compute q2: formula (A.6) in Berman and Eagleson (1985)
* Compute q3, q4 : Theorem (3.2) in Glaz (1989)
*---------------------------------------------------------*/
// compute the probability of term j = jmoy
q1 = 0.0;
for (j = 1; j <= jmoy; j++) {
jr = j;
q1 += Math.log (nr - jr + 1.0) - Math.log (jr);
}
q1 += jmoy*Math.log (d) + (nr - jmoy)*Math.log (q);
binMoy = Math.exp (q1);
q1 = binMoy;
jm1r = jmoy - m + 1;
if (((jmoy - m + 1) & 1) != 0)
signe = -1;
else
signe = 1;
q2 = signe*binMoy;
q3 = signe*binMoy*(2.0 - jm1r*jm1r + jm1r);
q4 = signe*binMoy*(jm1r + 1.0)*(jm1r + 2.0)*(6.0 + jm1r*jm1r -
5.0*jm1r);
// compute the probability of terms j > jmoy
if (((jmoy - m + 1) & 1) != 0)
signe = -1;
else
signe = 1;
jm1r = jmoy - m + 1;
bin = binMoy;
for (j = jmoy + 1; j <= n; j++) {
jr = j;
jm1r += 1.0;
signe = -signe;
bin = (bin*(nr - jr + 1.0)*d)/(jr*q);
if (bin < EPSILONSCAN)
break;
q1 += bin;
q2 += signe*bin;
q3 += signe*bin*(2.0 - jm1r*jm1r + jm1r);
q4 += signe*bin*(jm1r + 1.0)*(jm1r + 2.0)*(6.0 + jm1r*jm1r -
5.0*jm1r);
}
q1 = 1.0 - q1;
q3 /= 2.0;
q4 /= 12.0;
if (m == 3) {
// Problem with this formula; I do not get the same results as Glaz
q4 = ((nr*(nr - 1.0)*d*d*Math.pow (q, nr - 2.0))/8.0
+ nr*d*2.0*Math.pow (1.0 - 2.0*d, nr - 1.0))
- 4.0*Math.pow (1.0 - 2.0*d, nr);
if (d < 1.0/3.0) {
q4 += nr*d*2.0*Math.pow (1.0 - 3.0*d, nr - 1.0)
+ 4.0*Math.pow (1.0 - 3.0*d, nr);
}
}
// compute probability: Glaz, equations (3.2) and (3.3)
q3 = q1 - q2 - q3;
q4 = q3 - q4;
//when the approximation is bad, avoid overflow
temp = Math.log (q3) + (nr - m - 2.0)*Math.log (q4/q3);
if (temp >= 0.0)
return 0.0;
if (temp < -30.0)
return 1.0;
q4 = Math.exp (temp);
return 1.0 - q4;
}
//-----------------------------------------------------------------
private static double scanWNeff (int n, double d, int m) {
double q = 1.0 - d;
double temp;
double bin;
double sum;
int j;
/*--------------------------------------
* Anderson-Titterington: equation (4)
*--------------------------------------*/
// compute the probability of term j = m
sum = 0.0;
for (j = 1; j <= m; j++)
sum += Math.log ((double)(n - j + 1)) - Math.log ((double)j);
sum += m*Math.log (d) + (n - m)*Math.log (q);
bin = Math.exp (sum);
temp = (m/d - n - 1.0)*bin;
sum = bin;
// compute the probability of terms j > m
for (j = m + 1; j <= n; j++) {
bin *= (n - j + 1)*d/(j*q);
if (bin < EPSILONSCAN)
break;
sum += bin;
}
sum = 2.0*sum + temp;
return sum;
}
//----------------------------------------------------------------
private static double scanAsympt (int n, double d, int m) {
double kappa;
double temp;
double theta;
double sum;
/*--------------------------------------------------------------
* Anderson-Titterington: asymptotic formula after equation (4)
*--------------------------------------------------------------*/
theta = Math.sqrt (d/(1.0 - d));
temp = Math.sqrt ((double)n);
kappa = m/(d*temp) - temp;
temp = theta*kappa;
temp = temp*temp/2.0;
sum = 2.0*(1.0 - NormalDist.cdf01 (theta*kappa)) +
(kappa*theta*Math.exp (-temp))/(d*Math.sqrt (2.0*Math.PI));
return sum;
}
/**
* Return
* P[SN(d ) >= m], where SN(d ) is the scan
* statistic. It is defined as
*
*
*
* SN(d )= sup0 <= y <= 1-dη[y, y + d],
*
* where d is a constant in (0, 1),
*
* η[y, y + d] is the number of observations falling inside
* the interval [y, y + d], from a sample of N i.i.d. U(0, 1)
* random variables.
* The approximation returned by this function is generally good when
* it is close to 0, but is not very reliable when it exceeds, say, 0.4.
*
* Restrictions: N >= 2 and d <= 1/2.
*
* @param n sample size ( >= 2)
*
* @param d length of the test interval (∈(0, 1))
*
* @param m scan statistic
*
* @return the complementary distribution function of the statistic evaluated at m
*/
public static double scan (int n, double d, int m) {
double mu;
double prob;
if (n < 2)
throw new IllegalArgumentException ("Calling scan with n < 2");
if (d <= 0.0 || d >= 1.0)
throw new IllegalArgumentException ("Calling scan with "+
"d outside (0,1)");
if (m > n)
return 0.0;
if (m <= 1)
return 1.0;
if (m <= 2) {
if ((n - 1)*d >= 1.0)
return 1.0;
return 1.0 - Math.pow (1.0 - (n - 1)*d, (double)n);
}
if (d >= 0.5 && m <= (n + 1)/2.0)
return 1.0;
if (d > 0.5)
return -1.0; // Error
// util_Assert (d <= 0.5, "Calling fbar_Scan with d > 1/2");
mu = n*d; // mean of a binomial
if (m <= mu + d)
return 1.0;
if (mu <= 10.0)
return scanGlaz (n, d, m);
prob = scanAsympt (n, d, m);
if ((d >= 0.3 && n >= 50.0) || (n*d*d >= 250.0 && d < 0.3)) {
if (prob <= 0.4)
return prob;
}
prob = scanWNeff (n, d, m);
if (prob <= 0.4)
return prob;
prob = scanGlaz (n, d, m);
if (prob > 0.4 && prob <= 1.0)
return prob;
return 1.0;
}
}
Ads help maintain this website.