fix upper and lower incomplete gamma integral

src/shogun/mathematics/Statistics.cpp

@@ -555,6 +555,14 @@ float64_t CStatistics::gamma_incomplete_lower(float64_t s, float64_t z)
 	REQUIRE(s>0, "Given exponent (%f) must be greater than 0.\n", s);
 	REQUIRE(z>=0, "Given integral bound (%f) must be greater or equal to 0.\n", z);
 
+	/* save computation */
+	if (z<=0 || s<=0)
+		return 0;
+
+	/* is more numerically stable */
+	if (z>1 && z>s)
+		return 1-gamma_incomplete_upper(s, z);
+
 	// use series expansion
 	// inspiration from
 	// https://github.com/lh3/samtools/blob/master/bcftools/kfunc.c
@@ -571,26 +579,10 @@ float64_t CStatistics::gamma_incomplete_lower(float64_t s, float64_t z)
 	return CMath::exp(s * CMath::log(z) - z - CStatistics::lgamma(s + 1.) + CMath::log(sum));
 }
 
-float64_t CStatistics::lgamma_approx(float64_t z)
-{
-	float64_t x = 0;
-	x += 0.1659470187408462e-06 / (z+7);
-	x += 0.9934937113930748e-05 / (z+6);
-	x -= 0.1385710331296526     / (z+5);
-	x += 12.50734324009056      / (z+4);
-	x -= 176.6150291498386      / (z+3);
-	x += 771.3234287757674      / (z+2);
-	x -= 1259.139216722289      / (z+1);
-	x += 676.5203681218835      / z;
-	x += 0.9999999999995183;
-	return CMath::log(x) - 5.58106146679532777 - z + (z-0.5) * CMath::log(z+6.5);
-}
-
 float64_t CStatistics::gamma_incomplete_upper(float64_t s, float64_t z)
 {
 	REQUIRE(s>0, "Given exponent (%f) must be greater than 0.\n", s);
 	REQUIRE(z>=0, "Given integral bound (%f) must be greater or equal to 0.\n", z);
-
 	/* Taken inspiration from
 	 * https://github.com/lh3/samtools/blob/master/bcftools/kfunc.c
 	 * under MIT license
@@ -601,9 +593,23 @@ float64_t CStatistics::gamma_incomplete_upper(float64_t s, float64_t z)
 	 * (http://lib.stat.cmu.edu/apstat/245).
 	 */
 
+	/* save computation */
+	if (z<=0 || s<=0)
+		return 1;
+
+	/* is more numerically stable */
+	if (z<1 || z<s)
+		return 1-gamma_incomplete_lower(s, z);
+
+	/* save computation from underflows */
+	float64_t ax=s*CMath::log(z)-z-lgamma(s);
+	if (ax < -709.78271289338399)
+		return 0;
+
 	int32_t j;
 	float64_t C, D, f;
 	f = 1. + z - s;
+
 	C = f;
 	D = 0.;
 	// Modified Lentz's algorithm for computing continued fraction
@@ -625,7 +631,8 @@ float64_t CStatistics::gamma_incomplete_upper(float64_t s, float64_t z)
 		if (CMath::abs<float64_t>(d - 1.) < 1e-14)
 			break;
 	}
-	return CMath::exp(s * CMath::log(z) - z - CStatistics::lgamma_approx(s) - CMath::log(f));
+	float64_t result = CMath::exp(s * CMath::log(z) - z - CStatistics::lgamma(s) - CMath::log(f));
+	return result;
 }
 
 float64_t CStatistics::gamma_pdf(float64_t x, float64_t a, float64_t b)

src/shogun/mathematics/Statistics.h

@@ -230,31 +230,53 @@ class CStatistics: public CSGObject
 	 */
 	static float64_t gamma_pdf(float64_t x, float64_t a, float64_t b);
 
-	/** Compute the (regularized incomplete) Gamma function.
+	/** Compute the (lower) incomplete Gamma integral.
 	 *
-	 * Formulas are taken from Wiki, with additional input from Numerical
-	 * Recipes in C (for modified Lentz's algorithm). Uses a series expansion.
+	 * Given \f$p\f$, the function finds \f$x\f$ such that
+	 *
+	 * \f[
+	 * \text{incomplete\_gamma}(a,x)=\frac{1}{\Gamma(a)}}\int_0^x e^{-t} t^{a-1} dt.
+	 * \f]
+	 *
+	 *
+	 * In this implementation both arguments must be positive.
+	 * The integral is evaluated by either a power series or
+	 * continued fraction expansion, depending on the relative
+	 * values of \f$a\f$ and \f$x\f$.
 	 *
-	 * TODO latex
-	 * http://www.johndcook.com/blog/gamma_python/
+	 * Formulas are taken from Wiki, with additional input from Numerical
+	 * Recipes in C (for modified Lentz's algorithm for series expansion)
 	 *
-	 * Assumes a>0, x>=0
+	 * Assumes \f$a>0\f$, \f$x\geq 0\f$
 	 *
-	 * @param a TODO
-	 * @param x TODO
-	 * @return TODO
+	 * @param a Exponent parameter
+	 * @param x Lower integral bound
+	 * @return Value of incomplete Gamma integral
 	 */
 	static float64_t gamma_incomplete_lower(float64_t a, float64_t x);
 
-	/** Cumputes the (regularized upper) Gamma function.
+	/** Complemented incomplete Gamma integral (upper incomplete Gamma function)
+	 *
+	 * The function is defined by
+	 *
+	 * \f[
+	 * \text{incomplete\_gamma\_completed}(a,x)=1-\text{incomplete\_gamma}(a,x) =
+	 * \frac{1}{\Gamma (a)}\int_x^\infty e^{-t} t^{a-1} dt
+	 * \f]
+	 *
+	 * In this implementation both arguments must be positive.
+	 * The integral is evaluated by either a power series or
+	 * continued fraction expansion, depending on the relative
+	 * values of \f$a\f$ and \f$x\f$.
 	 *
 	 * Modified Lentz's algorithm for computing continued fraction.
 	 * See Numerical Recipes in C, 2nd edition, section 5.2.
 	 *
-	 * TODO latex
-	 * http://www.johndcook.com/blog/gamma_python/
+	 * Assumes \f$a>0\f$, \f$x\geq 0\f$
 	 *
-	 * Assumes a>0, x>=0
+	 * @param a Exponent parameter
+	 * @param x Lower integral bound
+	 * @return Value of incomplete Gamma integral
 	 */
 	static float64_t gamma_incomplete_upper(float64_t a, float64_t x);