|
0
|
1 #include <math.h>
|
|
|
2
|
|
|
3
|
|
|
4 /* Log gamma function
|
|
|
5 * \log{\Gamma(z)}
|
|
|
6 * AS245, 2nd algorithm, http://lib.stat.cmu.edu/apstat/245
|
|
|
7 */
|
|
|
8 double kf_lgamma(double z)
|
|
|
9 {
|
|
|
10 double x = 0;
|
|
|
11 x += 0.1659470187408462e-06 / (z+7);
|
|
|
12 x += 0.9934937113930748e-05 / (z+6);
|
|
|
13 x -= 0.1385710331296526 / (z+5);
|
|
|
14 x += 12.50734324009056 / (z+4);
|
|
|
15 x -= 176.6150291498386 / (z+3);
|
|
|
16 x += 771.3234287757674 / (z+2);
|
|
|
17 x -= 1259.139216722289 / (z+1);
|
|
|
18 x += 676.5203681218835 / z;
|
|
|
19 x += 0.9999999999995183;
|
|
|
20 return log(x) - 5.58106146679532777 - z + (z-0.5) * log(z+6.5);
|
|
|
21 }
|
|
|
22
|
|
|
23 /* complementary error function
|
|
|
24 * \frac{2}{\sqrt{\pi}} \int_x^{\infty} e^{-t^2} dt
|
|
|
25 * AS66, 2nd algorithm, http://lib.stat.cmu.edu/apstat/66
|
|
|
26 */
|
|
|
27 double kf_erfc(double x)
|
|
|
28 {
|
|
|
29 const double p0 = 220.2068679123761;
|
|
|
30 const double p1 = 221.2135961699311;
|
|
|
31 const double p2 = 112.0792914978709;
|
|
|
32 const double p3 = 33.912866078383;
|
|
|
33 const double p4 = 6.37396220353165;
|
|
|
34 const double p5 = .7003830644436881;
|
|
|
35 const double p6 = .03526249659989109;
|
|
|
36 const double q0 = 440.4137358247522;
|
|
|
37 const double q1 = 793.8265125199484;
|
|
|
38 const double q2 = 637.3336333788311;
|
|
|
39 const double q3 = 296.5642487796737;
|
|
|
40 const double q4 = 86.78073220294608;
|
|
|
41 const double q5 = 16.06417757920695;
|
|
|
42 const double q6 = 1.755667163182642;
|
|
|
43 const double q7 = .08838834764831844;
|
|
|
44 double expntl, z, p;
|
|
|
45 z = fabs(x) * M_SQRT2;
|
|
|
46 if (z > 37.) return x > 0.? 0. : 2.;
|
|
|
47 expntl = exp(z * z * - .5);
|
|
|
48 if (z < 10. / M_SQRT2) // for small z
|
|
|
49 p = expntl * ((((((p6 * z + p5) * z + p4) * z + p3) * z + p2) * z + p1) * z + p0)
|
|
|
50 / (((((((q7 * z + q6) * z + q5) * z + q4) * z + q3) * z + q2) * z + q1) * z + q0);
|
|
|
51 else p = expntl / 2.506628274631001 / (z + 1. / (z + 2. / (z + 3. / (z + 4. / (z + .65)))));
|
|
|
52 return x > 0.? 2. * p : 2. * (1. - p);
|
|
|
53 }
|
|
|
54
|
|
|
55 /* The following computes regularized incomplete gamma functions.
|
|
|
56 * Formulas are taken from Wiki, with additional input from Numerical
|
|
|
57 * Recipes in C (for modified Lentz's algorithm) and AS245
|
|
|
58 * (http://lib.stat.cmu.edu/apstat/245).
|
|
|
59 *
|
|
|
60 * A good online calculator is available at:
|
|
|
61 *
|
|
|
62 * http://www.danielsoper.com/statcalc/calc23.aspx
|
|
|
63 *
|
|
|
64 * It calculates upper incomplete gamma function, which equals
|
|
|
65 * kf_gammaq(s,z)*tgamma(s).
|
|
|
66 */
|
|
|
67
|
|
|
68 #define KF_GAMMA_EPS 1e-14
|
|
|
69 #define KF_TINY 1e-290
|
|
|
70
|
|
|
71 // regularized lower incomplete gamma function, by series expansion
|
|
|
72 static double _kf_gammap(double s, double z)
|
|
|
73 {
|
|
|
74 double sum, x;
|
|
|
75 int k;
|
|
|
76 for (k = 1, sum = x = 1.; k < 100; ++k) {
|
|
|
77 sum += (x *= z / (s + k));
|
|
|
78 if (x / sum < KF_GAMMA_EPS) break;
|
|
|
79 }
|
|
|
80 return exp(s * log(z) - z - kf_lgamma(s + 1.) + log(sum));
|
|
|
81 }
|
|
|
82 // regularized upper incomplete gamma function, by continued fraction
|
|
|
83 static double _kf_gammaq(double s, double z)
|
|
|
84 {
|
|
|
85 int j;
|
|
|
86 double C, D, f;
|
|
|
87 f = 1. + z - s; C = f; D = 0.;
|
|
|
88 // Modified Lentz's algorithm for computing continued fraction
|
|
|
89 // See Numerical Recipes in C, 2nd edition, section 5.2
|
|
|
90 for (j = 1; j < 100; ++j) {
|
|
|
91 double a = j * (s - j), b = (j<<1) + 1 + z - s, d;
|
|
|
92 D = b + a * D;
|
|
|
93 if (D < KF_TINY) D = KF_TINY;
|
|
|
94 C = b + a / C;
|
|
|
95 if (C < KF_TINY) C = KF_TINY;
|
|
|
96 D = 1. / D;
|
|
|
97 d = C * D;
|
|
|
98 f *= d;
|
|
|
99 if (fabs(d - 1.) < KF_GAMMA_EPS) break;
|
|
|
100 }
|
|
|
101 return exp(s * log(z) - z - kf_lgamma(s) - log(f));
|
|
|
102 }
|
|
|
103
|
|
|
104 double kf_gammap(double s, double z)
|
|
|
105 {
|
|
|
106 return z <= 1. || z < s? _kf_gammap(s, z) : 1. - _kf_gammaq(s, z);
|
|
|
107 }
|
|
|
108
|
|
|
109 double kf_gammaq(double s, double z)
|
|
|
110 {
|
|
|
111 return z <= 1. || z < s? 1. - _kf_gammap(s, z) : _kf_gammaq(s, z);
|
|
|
112 }
|
|
|
113
|
|
|
114 /* Regularized incomplete beta function. The method is taken from
|
|
|
115 * Numerical Recipe in C, 2nd edition, section 6.4. The following web
|
|
|
116 * page calculates the incomplete beta function, which equals
|
|
|
117 * kf_betai(a,b,x) * gamma(a) * gamma(b) / gamma(a+b):
|
|
|
118 *
|
|
|
119 * http://www.danielsoper.com/statcalc/calc36.aspx
|
|
|
120 */
|
|
|
121 static double kf_betai_aux(double a, double b, double x)
|
|
|
122 {
|
|
|
123 double C, D, f;
|
|
|
124 int j;
|
|
|
125 if (x == 0.) return 0.;
|
|
|
126 if (x == 1.) return 1.;
|
|
|
127 f = 1.; C = f; D = 0.;
|
|
|
128 // Modified Lentz's algorithm for computing continued fraction
|
|
|
129 for (j = 1; j < 200; ++j) {
|
|
|
130 double aa, d;
|
|
|
131 int m = j>>1;
|
|
|
132 aa = (j&1)? -(a + m) * (a + b + m) * x / ((a + 2*m) * (a + 2*m + 1))
|
|
|
133 : m * (b - m) * x / ((a + 2*m - 1) * (a + 2*m));
|
|
|
134 D = 1. + aa * D;
|
|
|
135 if (D < KF_TINY) D = KF_TINY;
|
|
|
136 C = 1. + aa / C;
|
|
|
137 if (C < KF_TINY) C = KF_TINY;
|
|
|
138 D = 1. / D;
|
|
|
139 d = C * D;
|
|
|
140 f *= d;
|
|
|
141 if (fabs(d - 1.) < KF_GAMMA_EPS) break;
|
|
|
142 }
|
|
|
143 return exp(kf_lgamma(a+b) - kf_lgamma(a) - kf_lgamma(b) + a * log(x) + b * log(1.-x)) / a / f;
|
|
|
144 }
|
|
|
145 double kf_betai(double a, double b, double x)
|
|
|
146 {
|
|
|
147 return x < (a + 1.) / (a + b + 2.)? kf_betai_aux(a, b, x) : 1. - kf_betai_aux(b, a, 1. - x);
|
|
|
148 }
|
|
|
149
|
|
|
150 #ifdef KF_MAIN
|
|
|
151 #include <stdio.h>
|
|
|
152 int main(int argc, char *argv[])
|
|
|
153 {
|
|
|
154 double x = 5.5, y = 3;
|
|
|
155 double a, b;
|
|
|
156 printf("erfc(%lg): %lg, %lg\n", x, erfc(x), kf_erfc(x));
|
|
|
157 printf("upper-gamma(%lg,%lg): %lg\n", x, y, kf_gammaq(y, x)*tgamma(y));
|
|
|
158 a = 2; b = 2; x = 0.5;
|
|
|
159 printf("incomplete-beta(%lg,%lg,%lg): %lg\n", a, b, x, kf_betai(a, b, x) / exp(kf_lgamma(a+b) - kf_lgamma(a) - kf_lgamma(b)));
|
|
|
160 return 0;
|
|
|
161 }
|
|
|
162 #endif
|
