/*************************************************************************\ * * 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 #include 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]; } /*=========================================================================*/