htestu-0.1.0.0: cbits/testu/src/fbar.c
/*************************************************************************\
*
* Package: ProbDist
* File: fbar.c
* Environment: ANSI C
*
* Copyright (c) 2002 Pierre L'Ecuyer, DIRO, Université de Montréal.
* e-mail: lecuyer@iro.umontreal.ca
* All rights reserved.
*
* Redistribution and use in source and binary forms, with or without
* modification, are permitted without a fee for private, research,
* academic, or other non-commercial purposes.
* Any use of this software in a commercial environment requires a
* written licence from the copyright owner.
*
* Any changes made to this package must be clearly identified as such.
*
* In scientific publications which used this software, a reference to it
* would be appreciated.
*
* Redistributions of source code must retain this copyright notice
* and the following disclaimer.
*
* THIS PACKAGE IS PROVIDED "AS IS" AND WITHOUT ANY EXPRESS OR
* IMPLIED WARRANTIES, INCLUDING, WITHOUT LIMITATION, THE IMPLIED
* WARRANTIES OF MERCHANTIBILITY AND FITNESS FOR A PARTICULAR PURPOSE.
*
\*************************************************************************/
#include "fbar.h"
#include "fdist.h"
#include "num.h"
#include "num2.h"
#include "util.h"
#include "gdef.h"
#include <math.h>
#include <float.h>
double fdist_belog(double);
extern const double fdist_XINF;
extern const double fdist_XBIG;
extern const double fdist_XBIGM;
/* EpsArray[j]: Epsilon required for j decimal degits of precision */
static const double EpsArray[] = {
0.5, 0.5E-1, 0.5E-2, 0.5E-3, 0.5E-4, 0.5E-5, 0.5E-6, 0.5E-7, 0.5E-8,
0.5E-9, 0.5E-10, 0.5E-11, 0.5E-12, 0.5E-13, 0.5E-14, 0.5E-15, 0.5E-16,
0.5E-17, 0.5E-18, 0.5E-19, 0.5E-20, 0.5E-21, 0.5E-22, 0.5E-23, 0.5E-24,
0.5E-25, 0.5E-26, 0.5E-27, 0.5E-28, 0.5E-29, 0.5E-30, 0.5E-31, 0.5E-32,
0.5E-33, 0.5E-34, 0.5E-35
};
/* Compute IMAX extra terms in the tails of discrete distributions */
static const long IMAX = 20;
/*=========================================================================*/
double fbar_Unif (double x)
{
if (x <= 0.0)
return 1.0;
if (x >= 1.0)
return 0.0;
return 1.0 - x;
}
/*=========================================================================*/
double fbar_Expon (double x)
{
if (x <= 0.0)
return 1.0;
if (x >= fdist_XBIGM)
return 0.0;
return exp (-x);
}
/*=========================================================================*/
double fbar_Weibull (double c, double x)
{
double temp;
util_Assert (c > 0.0, "fbar_Weibull: c <= 0");
if (x <= 0.0)
return 1.0;
if (x >= DBL_MAX_EXP * FLT_RADIX && c >= 1.0)
return 0.0;
temp = c*log(x);
if (temp >= DBL_MAX_EXP * num_Ln2)
return 0.0;
temp = exp(temp);
return (exp (-temp));
}
/*=========================================================================*/
double fbar_Logistic (double x)
{
if (x <= -fdist_XBIG) {
return 1.0;
}
if (x >= fdist_XBIG) {
return exp (-x);
}
return 1.0 / (1.0 + exp (x));
}
/*=========================================================================*/
double fbar_Pareto (double c, double x)
{
util_Assert (c > 0.0, "fbar_Pareto: c <= 0");
if (x <= 1.0)
return 1.0;
return (pow (x, -c));
}
/**************************************************************************/
double fbar_Normal1 (double x)
/*
* Returns P[X >= x] = 1 - F(x) where F is the normal distribution by
* computing the complementary distribution directly; it is thus more
* precise in the tail.
*/
{
static const double A[25] = {
6.10143081923200418E-1,
-4.34841272712577472E-1,
1.76351193643605501E-1,
-6.07107956092494149E-2,
1.77120689956941145E-2,
-4.32111938556729382E-3,
8.54216676887098679E-4,
-1.27155090609162743E-4,
1.12481672436711895E-5,
3.13063885421820973E-7,
-2.70988068537762022E-7,
3.07376227014076884E-8,
2.51562038481762294E-9,
-1.02892992132031913E-9,
2.99440521199499394E-11,
2.60517896872669363E-11,
-2.63483992417196939E-12,
-6.43404509890636443E-13,
1.12457401801663447E-13,
1.7281533389986098E-14,
-4.2641016949424E-15,
-5.4537197788E-16,
1.5869760776E-16,
2.08998378E-17,
-0.5900E-17
};
const double kk = 5.30330085889910643300; /* 3.75 Sqrt(2) */
double y, t;
int Neg;
if (x >= fdist_XBIG) {
return 0.0;
}
if (x <= -fdist_XBIG) {
return 1.0;
}
if (x >= 0.0)
Neg = 0;
else {
Neg = 1;
x = -x;
}
t = (x - kk) / (x + kk);
y = num2_EvalCheby (A, 24, t);
y = y * exp (-x * x / 2.0) / 2.0;
if (Neg == 1)
return (1.0 - y);
else
return (y);
}
/*=========================================================================*/
double fbar_Normal2 (double x)
{
static const double V[121] = {
1.2533141373155, 1.137490921203605, 1.037824575853727,
0.951527192071207, 0.8763644564536924, 0.8105337152790306,
0.7525711790634081, 0.7012808218544303, 0.6556795424187987,
0.61495459615093, 0.5784303460476312, 0.5455421356582171,
0.5158156382179634, 0.4888504415275737, 0.4643069280394423,
0.4418957328326002, 0.4213692292880546, 0.4025146181296722,
0.3851482907984348, 0.3691112106902635, 0.3542651113297938,
0.3404893532870847, 0.3276783146905521, 0.31573921586941,
0.3045902987101033, 0.2941592970402893, 0.284382146748493,
0.2752018941576065, 0.2665677689682238, 0.2584343943120386,
0.2507611114439651, 0.243511400615456, 0.2366523829135607,
0.230154390478801, 0.2239905946538289, 0.2181366833614714,
0.2125705804420318, 0.2072722008565011, 0.2022232366330547,
0.1974069692375194, 0.1928081047153158, 0.1884126285076003,
0.1842076773079702, 0.1801814257143918, 0.1763229857571027,
0.1726223176578506, 0.1690701504076941, 0.1656579109468773,
0.1623776608968675, 0.1592220399363674, 0.1561842150339759,
0.153257834853479, 0.1504369887362691, 0.1477161697413935,
0.145090241289131, 0.1425544070104023, 0.1401041834530503,
0.1377353753382303, 0.1354440530967635, 0.1332265324471292,
0.1310793558044918, 0.1289992753343376, 0.126983237485437,
0.1250283688553504, 0.1231319632579323, 0.1212914698765462,
0.119504482399253, 0.1177687290432979, 0.1160820633859823,
0.1144424559276431, 0.112847986320103, 0.1112968362007359,
0.1097872825783083, 0.1083176917221132, 0.1068865135106745,
0.1054922762005562, 0.1041335815795983, 0.1028091004723001,
0.1015175685681028, 0.1002577825460485, 0.09902859647173194,
0.09782891844465691, 0.09665770747608191, 0.09551397057921558,
0.09439676005522439, 0.09330517095996169, 0.09223833873763035,
0.09119543700877471, 0.09017567550106469, 0.08917829811230435,
0.08820258109597616, 0.08724783136042988, 0.08631338487354936,
0.08539860516539227, 0.08450288192189578, 0.08362562966329139,
0.08276628650136918, 0.08192431297018954, 0.08109919092525536,
0.08029042250654048, 0.07949752916111721, 0.07872005072144664,
0.07795754453568722, 0.07720958464664668, 0.07647576101624852,
0.07575567879261112, 0.07504895761704659, 0.07435523096847724,
0.07367414554294564, 0.07300536066605566, 0.07234854773633338,
0.07170338969763433, 0.07106958053885212, 0.07044682481930167,
0.06983483721825942, 0.06923334210724434, 0.06864207314371742,
0.06806077288496332, 0.0674891924209997, 0.06692709102543307,
0.06637423582325017
};
int j;
lebool negatif;
double t, u, z, h;
double r, r1, r2, r3, r4, r5, r6, r7, r8;
if (x >= fdist_XBIG) {
return 0.0;
}
if (x <= -fdist_XBIG) {
return 1.0;
}
if (x < 0.0) {
negatif = TRUE;
x = -x;
} else {
negatif = FALSE;
}
j = (int) (8.0 * x + 0.5);
if (j > 120)
j = 120;
z = 0.125 * j;
h = x - z;
r = V[j];
r1 = r * z - 1.0;
r2 = 0.5 * (r + z * r1);
r3 = (r1 + z * r2) / 3.0;
r4 = 0.25 * (r2 + z * r3);
r5 = 0.2 * (r3 + z * r4);
r6 = (r4 + z * r5) / 6.0;
r7 = (r5 + z * r6) / 7.0;
r8 = 0.125 * (r6 + z * r7);
t = r + h * (r1 + h * (r2 + h * (r3 + h * (r4 + h * (r5 + h * (r6 +
h * (r7 + h * r8)))))));
u = t * exp (-0.5 * x * x - 0.9189385332046727);
if (negatif)
return 1.0 - u;
else
return u;
}
/*=========================================================================*/
#ifdef HAVE_ERF
double fbar_Normal3 (double x)
{
return 0.5 * erfc (x * num_1Rac2);
}
#endif
/*=========================================================================*/
double fbar_BiNormal1 (double x, double y, double rho, int ndig)
{
return fdist_BiNormal1(-x, -y, rho, ndig);
}
/*=========================================================================*/
double fbar_BiNormal2 (double x, double y, double rho)
{
return fdist_BiNormal2 (-x, -y, rho);
}
/*=========================================================================*/
double fbar_LogNormal (double mu, double sigma, double x)
{
util_Assert (sigma > 0.0, "fbar_LogNormal: sigma <= 0");
if (x <= 0.0)
return 1.0;
return fbar_Normal1 ((log (x) - mu) / sigma);
}
/*=========================================================================*/
double fbar_JohnsonSB (double alpha, double beta, double a, double b,
double x)
{
util_Assert (beta > 0.0, "fbar_JohnsonSB: beta <= 0");
util_Assert (b > a, "fbar_JohnsonSB: b <= a");
if (x <= a)
return 1.0;
if (x >= b)
return 0.0;
return fbar_Normal1 (alpha + beta * log ((x - a) / (b - x)));
}
/*=========================================================================*/
double fbar_JohnsonSU (double alpha, double beta, double x)
{
const double XLIM = 1.0e10;
double r;
lebool negative = FALSE;
util_Assert (beta > 0.0, "fbar_JohnsonSU: beta <= 0");
if (x < 0.0) {
negative = TRUE;
x = -x;
}
/* compute r = x + sqrt (x * x + 1) */
if (x < XLIM)
r = x + sqrt (x * x + 1.0);
else
r = 2.0 * x;
if (negative)
r = 1.0 / r;
if (r > 0.0)
return fbar_Normal1 (alpha + beta * log (r));
else
return 1.0;
}
/*=========================================================================*/
static double F2AD[103]; /* Tables for the approximation of the */
static double CoAD[103]; /* Anderson-Darling distribution */
static void AndersonDarlingInit (void)
{
F2AD[0] = 0.0; F2AD[1] = 1.7315E-10;
F2AD[2] = 2.80781E-5; F2AD[3] = 1.40856E-3;
F2AD[4] = 9.58772E-3; F2AD[5] = 2.960552E-2;
F2AD[6] = 6.185146E-2; F2AD[7] = 1.0357152E-1;
F2AD[8] = 1.5127241E-1; F2AD[9] = 2.0190317E-1;
F2AD[10] = 2.5318023E-1; F2AD[11] = 3.0354278E-1;
F2AD[12] = 3.5200015E-1; F2AD[13] = 3.9797537E-1;
F2AD[14] = 4.4117692E-1; F2AD[15] = 4.8150305E-1;
F2AD[16] = 5.1897375E-1; F2AD[17] = 5.5368396E-1;
F2AD[18] = 5.8577199E-1; F2AD[19] = 6.1539864E-1;
F2AD[20] = 6.4273362E-1; F2AD[21] = 6.6794694E-1;
F2AD[22] = 6.9120359E-1; F2AD[23] = 7.126605E-1;
F2AD[24] = 7.3246483E-1; F2AD[25] = 7.507533E-1;
F2AD[26] = 7.6765207E-1; F2AD[27] = 7.8327703E-1;
F2AD[28] = 7.9773426E-1; F2AD[29] = 8.1112067E-1;
F2AD[30] = 8.2352466E-1; F2AD[31] = 8.3502676E-1;
F2AD[32] = 8.4570037E-1; F2AD[33] = 8.5561231E-1;
F2AD[34] = 8.6482346E-1; F2AD[35] = 8.7338931E-1;
F2AD[36] = 8.8136046E-1; F2AD[37] = 8.8878306E-1;
F2AD[38] = 8.9569925E-1; F2AD[39] = 9.0214757E-1;
F2AD[40] = 9.081653E-1; F2AD[41] = 9.1378043E-1;
F2AD[42] = 9.1902284E-1; F2AD[43] = 9.2392345E-1;
F2AD[44] = 9.2850516E-1; F2AD[45] = 9.3279084E-1;
F2AD[46] = 9.3680149E-1; F2AD[47] = 9.4055647E-1;
F2AD[48] = 9.440736E-1; F2AD[49] = 9.4736933E-1;
F2AD[50] = 9.5045883E-1; F2AD[51] = 9.5335611E-1;
F2AD[52] = 9.5607414E-1; F2AD[53] = 9.586249E-1;
F2AD[54] = 9.6101951E-1; F2AD[55] = 9.6326825E-1;
F2AD[56] = 9.6538067E-1; F2AD[57] = 9.6736563E-1;
F2AD[58] = 9.6923135E-1; F2AD[59] = 9.7098548E-1;
F2AD[60] = 9.7263514E-1; F2AD[61] = 9.7418694E-1;
F2AD[62] = 9.7564704E-1; F2AD[63] = 9.7702119E-1;
F2AD[64] = 9.7831473E-1; F2AD[65] = 9.7953267E-1;
F2AD[66] = 9.8067966E-1; F2AD[67] = 9.8176005E-1;
F2AD[68] = 9.827779E-1; F2AD[69] = 9.8373702E-1;
F2AD[70] = 9.8464096E-1; F2AD[71] = 9.8549304E-1;
F2AD[72] = 9.8629637E-1; F2AD[73] = 9.8705386E-1;
F2AD[74] = 9.8776824E-1; F2AD[75] = 9.8844206E-1;
F2AD[76] = 9.8907773E-1; F2AD[77] = 9.8967747E-1;
F2AD[78] = 9.9024341E-1; F2AD[79] = 9.9077752E-1;
F2AD[80] = 9.9128164E-1; F2AD[81] = 9.9175753E-1;
F2AD[82] = 9.9220682E-1; F2AD[83] = 9.9263105E-1;
F2AD[84] = 9.9303165E-1; F2AD[85] = 9.9340998E-1;
F2AD[86] = 9.9376733E-1; F2AD[87] = 9.9410488E-1;
F2AD[88] = 9.9442377E-1; F2AD[89] = 9.9472506E-1;
F2AD[90] = 9.9500974E-1; F2AD[91] = 9.9527876E-1;
F2AD[92] = 9.95533E-1; F2AD[93] = 9.9577329E-1;
F2AD[94] = 9.9600042E-1; F2AD[95] = 9.9621513E-1;
F2AD[96] = 9.964181E-1; F2AD[97] = 0.99661;
F2AD[98] = 9.9679145E-1; F2AD[99] = 9.9696303E-1;
F2AD[100] = 9.9712528E-1; F2AD[101] = 9.9727872E-1;
F2AD[102] = 9.9742384E-1;
CoAD[0] = 0.0;
CoAD[1] = 0.0; CoAD[2] = 0.0;
CoAD[3] = 0.0; CoAD[4] = 0.0;
CoAD[5] = -1.87E-3; CoAD[6] = 0.00898;
CoAD[7] = 0.0209; CoAD[8] = 0.03087;
CoAD[9] = 0.0377; CoAD[10] = 0.0414;
CoAD[11] = 0.04386; CoAD[12] = 0.043;
CoAD[13] = 0.0419; CoAD[14] = 0.0403;
CoAD[15] = 0.038; CoAD[16] = 3.54804E-2;
CoAD[17] = 0.032; CoAD[18] = 0.0293;
CoAD[19] = 2.61949E-2; CoAD[20] = 0.0228;
CoAD[21] = 0.0192; CoAD[22] = 1.59865E-2;
CoAD[23] = 0.0129; CoAD[24] = 0.0107;
CoAD[25] = 8.2464E-3; CoAD[26] = 0.00611;
CoAD[27] = 0.00363; CoAD[28] = 1.32272E-3;
CoAD[29] = -5.87E-4; CoAD[30] = -2.75E-3;
CoAD[31] = -3.95248E-3; CoAD[32] = -5.34E-3;
CoAD[33] = -6.892E-3; CoAD[34] = -8.10208E-3;
CoAD[35] = -8.93E-3; CoAD[36] = -9.552E-3;
CoAD[37] = -1.04605E-2; CoAD[38] = -0.0112;
CoAD[39] = -1.175E-2; CoAD[40] = -1.20216E-2;
CoAD[41] = -0.0124; CoAD[42] = -1.253E-2;
CoAD[43] = -1.27076E-2; CoAD[44] = -0.0129;
CoAD[45] = -1.267E-2; CoAD[46] = -1.22015E-2;
CoAD[47] = -0.0122; CoAD[48] = -1.186E-2;
CoAD[49] = -1.17218E-2; CoAD[50] = -0.0114;
CoAD[51] = -1.113E-2; CoAD[52] = -1.08459E-2;
CoAD[53] = -0.0104; CoAD[54] = -9.93E-3;
CoAD[55] = -9.5252E-3; CoAD[56] = -9.24E-3;
CoAD[57] = -9.16E-3; CoAD[58] = -8.8004E-3;
CoAD[59] = -8.63E-3; CoAD[60] = -8.336E-3;
CoAD[61] = -8.10512E-3; CoAD[62] = -7.94E-3;
CoAD[63] = -7.71E-3; CoAD[64] = -7.55064E-3;
CoAD[65] = -7.25E-3; CoAD[66] = -7.11E-3;
CoAD[67] = -6.834E-3; CoAD[68] = -0.0065;
CoAD[69] = -6.28E-3; CoAD[70] = -6.11008E-3;
CoAD[71] = -5.86E-3; CoAD[72] = -5.673E-3;
CoAD[73] = -5.35008E-3; CoAD[74] = -5.11E-3;
CoAD[75] = -4.786E-3; CoAD[76] = -4.59144E-3;
CoAD[77] = -4.38E-3; CoAD[78] = -4.15E-3;
CoAD[79] = -4.07696E-3; CoAD[80] = -3.93E-3;
CoAD[81] = -3.83E-3; CoAD[82] = -3.74656E-3;
CoAD[83] = -3.49E-3; CoAD[84] = -3.33E-3;
CoAD[85] = -3.20064E-3; CoAD[86] = -3.09E-3;
CoAD[87] = -2.93E-3; CoAD[88] = -2.78136E-3;
CoAD[89] = -2.72E-3; CoAD[90] = -2.66E-3;
CoAD[91] = -2.56208E-3; CoAD[92] = -2.43E-3;
CoAD[93] = -2.28E-3; CoAD[94] = -2.13536E-3;
CoAD[95] = -2.083E-3; CoAD[96] = -1.94E-3;
CoAD[97] = -1.82E-3; CoAD[98] = -1.77E-3;
CoAD[99] = -1.72E-3; CoAD[100] = -1.71104E-3;
CoAD[101] = -1.741E-3; CoAD[102] = -0.0016;
}
double fbar_AndersonDarling (long N, double X)
{
/* This function is not very precise for x < 0.05 */
const double h = 0.05; /* the step of the interpolation table */
static int ADFlag = 0;
double q;
double Res, Cor;
int i;
if (N == 1) {
if (X <= 0.38629436111989)
return 1.0;
if (X >= fdist_XBIGM)
return 0.0;
if (X < 6.0) {
q = 1.0 - 4.0 * exp(-X - 1.0);
return 1.0 - sqrt (q);
} else {
q = 4.0 * exp(-X - 1.0);
return 0.5*q*(1.0 + 0.25*q*(1.0 + 0.5*q*(1.0 + 0.125*q*(5.0 + 3.5*q))));
}
}
if (N <= 0) {
util_Warning (1, "fbar_AndersonDarling: N < 1");
return -1.0;
}
if (X > 10.0)
/* Sinclair-Spurr upper tail approximation (3.5) */
return 1.732 * exp(-X) / sqrt(num_Pi * X);
if (X > 5.0) {
/* asymptotic X: our empirical fit */
Res = exp (-0.56 - 1.06 * X);
q = exp (-1.03 - 1.06 * X); /* Empirical correction in 1/N */
return Res + q / N;
}
if (X <= 0.2)
return 1.0 - fdist_AndersonDarling (N, X);
if (ADFlag == 0) {
AndersonDarlingInit ();
ADFlag = 1;
}
i = 1 + (int) (X / h);
q = X / h - i;
/* Newton backwards quadratic interpolation */
Res = (F2AD[i - 2] - 2.0 * F2AD[i - 1] + F2AD[i]) * q * (q + 1.0) / 2.0
+ (F2AD[i] - F2AD[i - 1]) * q + F2AD[i];
/* Empirical correction in 1/N */
Cor = (CoAD[i] * (q + 1.0) - CoAD[i - 1] * q) / N;
Res = 1.0 - Res - Cor;
if (Res >= 1.0)
return 1.0;
if (Res <= 0.0)
return 0.0;
return Res;
}
/*=========================================================================*/
double fbar_ChiSquare1 (long N, double x)
/*
* Returns an approximation of the complementary Chi square cdf (N degrees
* of freedom). Similar to p:116 of W.J.Kennedy Jr and J.E.Gentle.
* Statistical computing, Dekker, New York, 1980. More precise in the
* tail than simply returning 1 - fdist_ChiSquare.
*/
{
const double XBIG_CHI = 2000.0;
const double tiers = 0.33333333333333333;
const double pt2 = 0.22222222222222222;
const double moinshuit = -8.3;
const double gam = 0.8862269254527579825931;
double H, E, DemiX, Terme, Sommation, Y;
long i;
util_Assert (N > 0, "Calling fbar_ChiSquare1 with N < 1");
if (x <= 0.0)
return 1.0;
if (N >= 150) {
if (x >= N * fdist_XBIG)
return 0.0;
} else {
if (x >= XBIG_CHI)
return 0.0;
}
if (N > 1000) {
if (x < 2.0)
return 1.0;
x = (pow ((x / N), tiers) - (1.0 - pt2 / N)) / sqrt (pt2 / N);
if (x > 35.0)
return 0.0;
if (x <= moinshuit)
return 1.0;
return fbar_Normal1 (x);
}
DemiX = x / 2.0;
if (!(N & 1)) { /* even N */
Terme = exp (-DemiX);
Sommation = Terme;
for (i = 1; i < N / 2; i++) {
Terme = Terme * DemiX / i;
Sommation += Terme;
}
Y = Sommation;
} else {
H = 2.0 * fbar_Normal1 (sqrt (x));
if (N == 1)
return H;
E = exp (-DemiX);
Terme = sqrt (DemiX) * E / gam;
for (i = 3; i < N; i += 2) {
H += Terme;
Terme = Terme * DemiX * 2.0 / i;
}
Y = H + Terme;
}
if (Y > 1.0)
return 1.0;
else
return Y;
}
/*=========================================================================*/
double fbar_ChiSquare2 (long n, int d, double x)
{
util_Assert (n > 0, "fbar_ChiSquare2: n <= 0");
if (x <= 0.0)
return 1.0;
return fbar_Gamma (n / 2.0, d, x / 2.0);
}
/*=========================================================================*/
double fbar_Gamma (double alpha, int d, double x)
{
const double aLIM = 1.0E5;
const double RENORM = 1.0E100;
const double EPS = EpsArray[d];
double V[6];
double v, res, A, B, R, term, dif;
int i;
util_Assert (alpha > 0.0, "fbar_Gamma: a <= 0");
util_Assert (d > 0, "fbar_Gamma: d <= 0");
util_Assert (d < 16, "fbar_Gamma: d > 15");
if (x <= 0.0)
return 1.0;
if (1.0 == alpha)
return fbar_Expon (x);
if (alpha >= 70.0) {
if (x >= alpha * fdist_XBIG)
return 0.0;
} else {
if (x >= fdist_XBIGM)
return 0.0;
}
if (alpha >= aLIM) {
double d2 = x + 1.0/3.0 - alpha - 0.02/alpha;
double S = alpha - 1.0/2.0;
double z = d2 * sqrt((1 + fdist_belog(S/x))/x);
return fbar_Normal1 (z);
}
if (x <= 1.0 || x < alpha)
return 1.0 - fdist_Gamma (alpha, d, x);
v = exp (alpha * log (x) - x - num2_LnGamma (alpha));
A = 1.0 - alpha;
B = A + x + 1.0;
term = 0.0;
V[0] = 1.0;
V[1] = x;
V[2] = x + 1.0;
V[3] = x * B;
res = V[2] / V[3];
do {
A += 1.0;
B += 2.0;
term += 1.0;
V[4] = B * V[2] - A * term * V[0];
V[5] = B * V[3] - A * term * V[1];
if (V[5] != 0.0) {
R = V[4] / V[5];
dif = fabs (res - R);
if (dif <= EPS * R)
return (v * res);
res = R;
}
for (i = 0; i < 4; i++)
V[i] = V[i + 2];
if (fabs (V[4]) >= RENORM) {
for (i = 0; i < 4; i++)
V[i] /= RENORM;
}
} while (1);
/* to eliminate a warning from the compiler; never reached */
return 0.0;
}
/*=========================================================================*/
static double KSPlusbarAsymp (long n, double x)
{
/* Compute the probability of the 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*exp(-z);
if (v >= 1.0)
return 1.0;
return v;
}
/*-------------------------------------------------------------------------*/
static double KSPlusbarUpper (long n, double x)
{
/* Compute the probability of the KSPlus distribution in the upper
tail using Smirnov's stable formula */
const double EPSILON = 1.0E-10;
double q;
double Sum = 0.0;
double term;
double t;
double LogCom;
double LOGJMAX;
int j;
int jmax = (int)(n - n*x);
/* We must avoid log(0) for j = jmax and q ~ 1.0 */
if ((1.0 - x - (double)jmax / n) <= 0.0)
jmax--;
j = jmax/2;
LogCom = num2_LnFactorial((int)n) - num2_LnFactorial(j) -
num2_LnFactorial((int)(n-j));
LOGJMAX = LogCom;
while (j > 0) {
q = (double)j / n + x;
term = LogCom + (j - 1)*log (q) + (n - j)*num2_log1p (-q);
t = exp (term);
Sum += t;
LogCom += log ((double)j / (n - j + 1));
if (t <= Sum*EPSILON)
break;
j--;
}
j = jmax/2;
LogCom = LOGJMAX + log ((double)(n - j)/(j + 1));
j++;
while (j <= jmax) {
q = (double)j / n + x;
term = LogCom + (j - 1)*log(q) + (n - j)*num2_log1p(-q);
t = exp (term);
Sum += t;
LogCom += log ((double)(n - j)/(j + 1));
if (t <= Sum*EPSILON)
break;
j++;
}
Sum *= x;
/* add the term j = 0 */
Sum += exp (n*num2_log1p (-x));
return Sum;
}
/*=========================================================================*/
double fbar_KSPlus (long N, double x)
{
const double NxParam = 6.5; /* frontier: alternating series */
const long NParam = 4000; /* frontier: non-alternating series */
const long NAsymp = 200000; /* frontier: asymptotic */
util_Assert (N > 0, "Calling fbar_KSPlus with N < 1");
if (x <= 0.0)
return 1.0;
if ((x >= 1.0) || (N*x*x >= 370.0))
return 0.0;
if (N == 1)
return 1.0 - x;
if (N * x <= NxParam)
return 1.0 - fdist_KSPlus (N, x);
if (N >= NAsymp)
return KSPlusbarAsymp (N, x);
if ((N <= NParam) || (N*x*x > 1.0))
return KSPlusbarUpper(N, x);
/* return (1.0 - 2.0*x/3.0)*exp(-2.0*N*x*x); */
return KSPlusbarAsymp (N, x);
}
/*=========================================================================*/
static double KSSpecial (long n, double x)
{
#define NLIM 20
if ((n * x * x >= 370.0) || (x >= 1.0))
return 0.0;
if (x <= 0.5 / n)
return 1.0;
if (n == 1)
return 2.0 - 2.0 * x;
if (x <= 1.0 / n) {
double w;
double t = 2.0 * x - 1.0 / n;
if (n <= NLIM) {
w = num2_Factorial ((int) n);
return 1.0 - w * pow (t, (double) n);
}
w = num2_LnFactorial ((int) n) + n * log (t);
return 1.0 - exp (w);
}
if (x >= 1.0 - 1.0 / n) {
return 2.0 * pow (1.0 - x, (double) n);
}
return -1.0;
}
#undef NLIM
/*-------------------------------------------------------------------------*/
double fbar_KS1 (long n, double x)
{
double v = KSSpecial(n, x);
if (v >= 0.0)
return v;
if (n <= 400) {
if (n*x*x < 4.0)
return 1.0 - fdist_KS1(n, x);
else
return 2.0 * KSPlusbarUpper(n, x);
}
if (n*x*x >= 2.2) {
if (n <= 200000)
return 2.0 * KSPlusbarUpper(n, x);
return 2.0*KSPlusbarAsymp (n, x);
}
return 1.0 - fdist_KS1(n, x);
}
/*=========================================================================*/
double fbar_CramerMises (long N, double x)
{
return 1.0 - fdist_CramerMises (N, x);
}
double fbar_WatsonG (long N, double x)
{
return 1.0 - fdist_WatsonG (N, x);
}
/*=========================================================================*/
double fbar_WatsonU (long N, double x)
{
/*
* Only the asymptotic form has been implemented. In the trivial case
* N = 1, we simply return 0.5
*/
const double xSepare = 0.15;
if (x <= 0.0)
return 1.0;
if (x >= fdist_XBIG)
return 0.0;
if (N == 1) /* N = 1, degenerate case */
return 0.5;
if (x > xSepare) {
/* this series converges rapidly for x > 0.15 */
const int JMAX = 10;
int j;
double signe;
double v;
double terme;
double somme;
v = exp (-(x * 2.0 * num_Pi * num_Pi));
signe = 1.0;
somme = 0.0;
j = 1;
do {
terme = pow (v, (double) j * j);
somme += signe * terme;
signe = -signe;
++j;
} while (!(terme < DBL_EPSILON || j > JMAX));
util_Warning (j > JMAX, "fbar_WatsonU: sum1 has not converged");
v = 2.0 * somme;
if (v <= 0.0)
return 0.0;
return v;
}
return 1.0 - fdist_WatsonU (N, x);
}
/*=========================================================================*/
/******************************\
*
* DISCRETE DISTRIBUTIONS
*
\******************************/
/*=========================================================================*/
static const double epsilonScan = 1.0E-7;
static double ScanGlaz (long N, double d, long m)
{
long 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 = (long) ((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 += log (Nr - jr + 1.0) - log (jr);
}
Q1 += jmoy * log (d) + (Nr - jmoy) * log (q);
BinMoy = exp (Q1);
Q1 = BinMoy;
jm1r = jmoy - m + 1;
if ((jmoy - m + 1) & 1)
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)
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 * pow (q, Nr - 2.0)) / 8.0
+ Nr * d * 2.0 * pow (1.0 - 2.0 * d, Nr - 1.0))
- 4.0 * pow (1.0 - 2.0 * d, Nr);
if (d < 1.0 / 3.0) {
Q4 += Nr * d * 2.0 * pow (1.0 - 3.0 * d, Nr - 1.0)
+ 4.0 * 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 = log (Q3) + (Nr - m - 2.0) * log (Q4 / Q3);
if (temp >= 0.0)
return 0.0;
if (temp < (-30.0))
return 1.0;
Q4 = exp (temp);
return 1.0 - Q4;
}
/*----------------------------------------------------------------------*/
static double ScanWNeff (long N, double d, long m)
{
double q = 1.0 - d;
double temp;
double Bin;
double Sum;
long j;
/*--------------------------------------*/
/* Anderson-Titterington: equation (4) */
/*--------------------------------------*/
/* compute the probability of term j = m */
Sum = 0.0;
for (j = 1; j <= m; j++) {
Sum += log ((double) (N - j + 1)) - log ((double) j);
}
Sum += m * log (d) + (N - m) * log (q);
Bin = 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;
}
/*----------------------------------------------------------------------*/
static double ScanAsympt (long N, double d, long m)
{
double Kappa;
double temp;
double Theta;
double Sum;
/*--------------------------------------------------------------*/
/* Anderson-Titterington: asymptotic formula after equation (4) */
/*--------------------------------------------------------------*/
Theta = sqrt (d / (1.0 - d));
temp = sqrt ((double) N);
Kappa = m / (d * temp) - temp;
temp = Theta * Kappa;
temp = temp * temp / 2.0;
Sum = 2.0 * fbar_Normal1 (Theta * Kappa) +
(Kappa * Theta * exp (-temp)) / (d * sqrt (2.0 * num_Pi));
return Sum;
}
/*----------------------------------------------------------------------*/
double fbar_Scan (long N, double d, long m)
{
double mu;
double prob;
util_Assert (N >= 2, "Calling fbar_Scan with N < 2");
util_Assert (d > 0.0 && d < 1.0,
"Calling fbar_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 - 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;
}
/*=========================================================================*/
double fbar_Geometric (double p, long n)
{
util_Assert (p >= 0.0 && p <= 1.0, "fbar_Geometric: p not in [0, 1]");
if (n <= 0)
return 1.0;
if (p >= 1.0) /* In fact, p == 1 */
return 0.0;
if (p <= 0.0) /* In fact, p == 0 */
return 1.0;
return pow (1.0 - p, (double) n);
}
/*=========================================================================*/
double fbar_Poisson1 (double lam, long s)
{
const double lamlim = 150.0;
long i;
double term, sum;
util_Assert (lam >= 0.0, "fbar_Poisson1: lambda < 0");
if (s <= 0)
return 1.0;
/* If lam > lamlim, we use the Chi2 distribution according to the exact
relation, with 2s + 2 degrees of freedom
fdist_Poisson (lam, s) = 1 - fdist_ChiSquare (2s + 2, 2*lam)
which also equals 1 - fdist_Gamma (s + 1, lam) */
if (lam > lamlim)
return fdist_Gamma ((double) s, 15, lam);
if (s <= lam)
return 1.0 - fdist_Poisson1 (lam, s - 1);
/* Sum at least IMAX prob. terms from i = s to i = oo */
sum = term = fmass_PoissonTerm1 (lam, s);
i = s + 1;
while (term > fmass_Epsilon || i <= s + IMAX) {
term *= lam / i;
sum += term;
i++;
}
return sum;
}
/*=========================================================================*/
double fbar_Poisson2 (fmass_INFO W, long s)
/*
* fbar_Poisson (lam, s) = 1 - fdist_Poisson (lam, s - 1)
*/
{
double lam;
util_Assert (W != NULL, "fbar_Poisson2: fmass_INFO is NULL pointer");
lam = W->paramR[0];
if (s <= 0)
return 1.0;
/* For large lam, we use the Chi2 distribution according to the exact
relation, with 2s + 2 degrees of freedom
fdist_Poisson (lam, s) = 1 - fdist_ChiSquare (2s + 2, 2*lam)
fdist_Poisson (lam, s) = 1 - fdist_Gamma (s + 1, lam)
*/
if (W->cdf == NULL)
return fdist_Gamma ((double) s, 15, lam);
if (s > W->smax)
return fbar_Poisson1 (lam, s);
if (s < W->smin)
return 1.0;
if (s > W->smed)
/* We keep the complementary distribution in the upper part of cdf */
return W->cdf[s - W->smin];
else
return 1.0 - W->cdf[s - 1 - W->smin];
}
/*=========================================================================*/
double fbar_Binomial2 (fmass_INFO W, long s)
{
double p;
long n;
util_Assert (W != NULL, "fbar_Binomial2: fmass_INFO is NULL pointer");
n = W->paramI[0];
p = W->paramR[0];
util_Assert (p >= 0.0 && p <= 1.0, "fbar_Binomial2: p not in [0, 1]");
if (0 == n)
return 1.0;
if (s < 1)
return 1.0;
if (s > n)
return 0.0;
if (p == 0.0)
return 0.0;
if (p == 1.0)
return 1.0;
if (W->cdf != NULL) {
if (s >= W->smax) {
/* Add IMAX dominant terms to get a few decimals in the tail */
const double q = 1.0 - p;
double z, sum, term;
long i;
sum = term = fmass_BinomialTerm3 (n, p, s);
if (fabs (q) > 0.0) {
z = p / q;
} else {
z = 0.0;
util_Warning (1, "fbar_Binomial2: p / q = infinite");
}
i = s;
while (i < n && i < s + IMAX) {
term = term * z * (n - i) / (i + 1);
sum += term;
i++;
}
return sum;
/* return fdist_Beta (s, n - s + 1, 10, p); */
}
if (s <= W->smin)
return 1.0;
if (s > W->smed)
/* We keep the complementary distribution in the upper part of cdf */
return W->cdf[s - W->smin];
else
return 1.0 - W->cdf[s - 1 - W->smin];
} else {
return 1.0 - fdist_Binomial1 (n, p, s - 1);
}
}
/*=========================================================================*/
double fbar_NegaBin2 (fmass_INFO W, long s)
{
double p;
long n;
util_Assert (W != NULL, "fbar_NegaBin2: fmass_INFO is NULL pointer");
n = W->paramI[0];
p = W->paramR[0];
util_Assert (p >= 0.0 && p <= 1.0, "fbar_NegaBin2: p not in [0, 1]");
if (s < 1)
return 1.0;
if (p >= 1.0) /* In fact, p == 1 */
return 0.0;
if (p <= 0.0) /* In fact, p == 0 */
return 1.0;
if (W->cdf == NULL)
return fdist_Binomial1 (s - 1 + n, p, n - 1);
if (s >= W->smax)
return fdist_Binomial1 (s - 1 + n, p, n - 1);
if (s <= W->smin)
return 1.0;
if (s > W->smed)
/* We keep the complementary distribution in the upper part of cdf */
return W->cdf[s - W->smin];
else
return 1.0 - W->cdf[s - 1 - W->smin];
}
/*=========================================================================*/