X-Git-Url: https://git.donarmstrong.com/?a=blobdiff_plain;f=bcftools%2Fkfunc.c;h=a637b6ca70afd768e736117b0dbbc377d36d2318;hb=0b3418cf166ce4a58cedf0d9a2df5ec3dd4cc5fa;hp=e9d7318ef482bcbce1020ded40e3095e21d417ca;hpb=6d8992db11070d00f712916c9f7dee4e49ff6efa;p=samtools.git

diff --git a/bcftools/kfunc.c b/bcftools/kfunc.c
index e9d7318..a637b6c 100644
--- a/bcftools/kfunc.c
+++ b/bcftools/kfunc.c
@@ -52,80 +52,111 @@ double kf_erfc(double x)
 	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