return x > 0.? 2. * p : 2. * (1. - p);
}
-/* Regularized (incomplete lower) gamma function
- * \frac{\gamma(p,x)}{\Gamma(p)}=\frac{1}{\Gamma(p)} \int_0^x t^{p-1}e^{-t} dt
- * AS245, http://lib.stat.cmu.edu/apstat/245
+/* The following computes regularized incomplete gamma functions.
+ * Formulas are taken from Wiki, with additional input from Numerical
+ * Recipes in C (for modified Lentz's algorithm) and AS245
+ * (http://lib.stat.cmu.edu/apstat/245).
+ *
+ * A good online calculator is available at:
+ *
+ * http://www.danielsoper.com/statcalc/calc23.aspx
+ *
+ * It calculates upper incomplete gamma function, which equals
+ * kf_gammaq(s,z)*tgamma(s).
*/
-double kf_gammap(double p, double x)
+
+#define KF_GAMMA_EPS 1e-14
+#define KF_TINY 1e-290
+
+// regularized lower incomplete gamma function, by series expansion
+static double _kf_gammap(double s, double z)
+{
+ double sum, x;
+ int k;
+ for (k = 1, sum = x = 1.; k < 100; ++k) {
+ sum += (x *= z / (s + k));
+ if (x / sum < KF_GAMMA_EPS) break;
+ }
+ return exp(s * log(z) - z - kf_lgamma(s + 1.) + log(sum));
+}
+// regularized upper incomplete gamma function, by continued fraction
+static double _kf_gammaq(double s, double z)
{
- double ret_val;
- double a, b, c, an, rn, pn1, pn2, pn3, pn4, pn5, pn6, arg;
+ int j;
+ double C, D, f;
+ f = 1. + z - s; C = f; D = 0.;
+ // Modified Lentz's algorithm for computing continued fraction
+ // See Numerical Recipes in C, 2nd edition, section 5.2
+ for (j = 1; j < 100; ++j) {
+ double a = j * (s - j), b = (j<<1) + 1 + z - s, d;
+ D = b + a * D;
+ if (D < KF_TINY) D = KF_TINY;
+ C = b + a / C;
+ if (C < KF_TINY) C = KF_TINY;
+ D = 1. / D;
+ d = C * D;
+ f *= d;
+ if (fabs(d - 1.) < KF_GAMMA_EPS) break;
+ }
+ return exp(s * log(z) - z - kf_lgamma(s) - log(f));
+}
- if (x == 0.) return 0.;
- // The following line is not thoroughly tested, so it is commented out.
- if (p > 1e3) return .5 * kf_erfc(-M_SQRT1_2 * sqrt(p) * 3. * (pow(x / p, 1./3.) + 1. / (p * 9.) - 1.));
- if (x > 1e8) return 1.;
- if (x <= 1. || x < p) { // series expansion
- c = 1.;
- arg = p * log(x) - x - kf_lgamma(p + 1.);
- ret_val = 1.;
- a = p;
- while (c > 1e-14) {
- a += 1.;
- c = c * x / a;
- ret_val += c;
- }
- arg += log(ret_val);
- ret_val = 0.;
- if (arg >= -88.) ret_val = exp(arg);
- } else { // continued expansion
- arg = p * log(x) - x - kf_lgamma(p);
- a = 1. - p;
- b = a + x + 1.;
- c = 0.;
- pn1 = 1.;
- pn2 = x;
- pn3 = x + 1.;
- pn4 = x * b;
- ret_val = pn3 / pn4;
- while (1) {
- a += 1.; b += 2.; c += 1.;
- an = a * c;
- pn5 = b * pn3 - an * pn1;
- pn6 = b * pn4 - an * pn2;
- if (fabs(pn6) > 0.) {
- rn = pn5 / pn6;
- if (fabs(ret_val - rn) <= fmin(1e-14, rn * 1e-14)) break;
- ret_val = rn;
- }
- pn1 = pn3; pn2 = pn4; pn3 = pn5; pn4 = pn6;
- if (fabs(pn5) >= 1e37)
- pn1 /= 1e37, pn2 /= 1e37, pn3 /= 1e37, pn4 /= 1e37;
- }
- arg += log(ret_val);
- ret_val = 1.;
- if (arg >= -88.) ret_val = 1. - exp(arg);
- }
- return ret_val;
+double kf_gammap(double s, double z)
+{
+ return z <= 1. || z < s? _kf_gammap(s, z) : 1. - _kf_gammaq(s, z);
}
-/* Numerical Recipe separates series expansion and continued
- * expansion. This may potentially reduce underflow for some
- * combinations of p and x. Nonetheless, the precision here is good
- * enough for me. I will not spend more time for now.
+double kf_gammaq(double s, double z)
+{
+ return z <= 1. || z < s? 1. - _kf_gammap(s, z) : _kf_gammaq(s, z);
+}
+
+/* Regularized incomplete beta function. The method is taken from
+ * Numerical Recipe in C, 2nd edition, section 6.4. The following web
+ * page calculates the incomplete beta function, which equals
+ * kf_betai(a,b,x) * gamma(a) * gamma(b) / gamma(a+b):
+ *
+ * http://www.danielsoper.com/statcalc/calc36.aspx
*/
-double kf_gammaq(double p, double x)
+static double kf_betai_aux(double a, double b, double x)
+{
+ double C, D, f;
+ int j;
+ if (x == 0.) return 0.;
+ if (x == 1.) return 1.;
+ f = 1.; C = f; D = 0.;
+ // Modified Lentz's algorithm for computing continued fraction
+ for (j = 1; j < 200; ++j) {
+ double aa, d;
+ int m = j>>1;
+ aa = (j&1)? -(a + m) * (a + b + m) * x / ((a + 2*m) * (a + 2*m + 1))
+ : m * (b - m) * x / ((a + 2*m - 1) * (a + 2*m));
+ D = 1. + aa * D;
+ if (D < KF_TINY) D = KF_TINY;
+ C = 1. + aa / C;
+ if (C < KF_TINY) C = KF_TINY;
+ D = 1. / D;
+ d = C * D;
+ f *= d;
+ if (fabs(d - 1.) < KF_GAMMA_EPS) break;
+ }
+ return exp(kf_lgamma(a+b) - kf_lgamma(a) - kf_lgamma(b) + a * log(x) + b * log(1.-x)) / a / f;
+}
+double kf_betai(double a, double b, double x)
{
- return 1. - kf_gammap(p, x);
+ return x < (a + 1.) / (a + b + 2.)? kf_betai_aux(a, b, x) : 1. - kf_betai_aux(b, a, 1. - x);
}
#ifdef KF_MAIN
#include <stdio.h>
int main(int argc, char *argv[])
{
- double x = 10, y = 2.5;
+ double x = 5.5, y = 3;
+ double a, b;
printf("erfc(%lg): %lg, %lg\n", x, erfc(x), kf_erfc(x));
- printf("lower-gamma(%lg,%lg): %lg\n", x, y, (1.0-kf_gammap(y, x))*tgamma(y));
+ printf("upper-gamma(%lg,%lg): %lg\n", x, y, kf_gammaq(y, x)*tgamma(y));
+ a = 2; b = 2; x = 0.5;
+ 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)));
return 0;
}
#endif
projects / samtools.git / blobdiff