diff --git a/src/shogun/mathematics/Statistics.cpp b/src/shogun/mathematics/Statistics.cpp
index 65f05a6e583..b712bd75fcb 100644
--- a/src/shogun/mathematics/Statistics.cpp
+++ b/src/shogun/mathematics/Statistics.cpp
@@ -661,6 +661,12 @@ float64_t CStatistics::inverse_gamma_cdf(float64_t p, float64_t a,
 	return inverse_incomplete_gamma_completed(a, 1-p)*b;
 }
 
+float64_t CStatistics::incomplete_gamma(float64_t a, float64_t x)
+{
+	SG_SERROR("NOT IMPLEMENTED");
+	return 0;
+}
+
 //float64_t CStatistics::incomplete_beta(float64_t a, float64_t b, float64_t x)
 //{
 //	SG_SERROR("NOT IMPLEMENTED");
diff --git a/src/shogun/mathematics/Statistics.h b/src/shogun/mathematics/Statistics.h
index a2d65878da8..14e73aa1d7a 100644
--- a/src/shogun/mathematics/Statistics.h
+++ b/src/shogun/mathematics/Statistics.h
@@ -214,6 +214,24 @@ class CStatistics: public CSGObject
 		return ::tgamma((double) x);
 	}
 
+	/** Incomplete gamma integral
+	 *
+	 * 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$.
+	 *
+	 * Taken from ALGLIB under gpl2+
+	 */
+	static float64_t incomplete_gamma(float64_t a, float64_t x);
+
 	/** Evaluates the CDF of the gamma distribution with given parameters \f$a, b\f$
 	 * at \f$x\f$.
 	 *