diff --git a/src/Numerics/SpecialFunctions/Gamma.cs b/src/Numerics/SpecialFunctions/Gamma.cs
index c5c21d46..ff2ace6f 100644
--- a/src/Numerics/SpecialFunctions/Gamma.cs
+++ b/src/Numerics/SpecialFunctions/Gamma.cs
@@ -70,7 +70,7 @@ namespace MathNet.Numerics
};
///
- /// Computes the logarithm of the Gamma function.
+ /// Computes the logarithm of the Gamma function.
///
/// The argument of the gamma function.
/// The logarithm of the gamma function.
@@ -112,7 +112,7 @@ namespace MathNet.Numerics
}
///
- /// Computes the Gamma function.
+ /// Computes the Gamma function.
///
/// The argument of the gamma function.
/// The logarithm of the gamma function.
@@ -151,7 +151,7 @@ namespace MathNet.Numerics
return s * Constants.TwoSqrtEOverPi * Math.Pow((z - 0.5 + GammaR) / Math.E, z - 0.5);
}
}
-
+
///
/// Returns the upper incomplete regularized gamma function
/// Q(a,x) = 1/Gamma(a) * int(exp(-t)t^(a-1),t=0..x) for real a > 0, x > 0.
@@ -164,7 +164,7 @@ namespace MathNet.Numerics
const double epsilon = 0.000000000000001;
const double big = 4503599627370496.0;
const double bigInv = 2.22044604925031308085e-16;
-
+
if (x <= 0d || a <= 0d)
{
return 1d;
@@ -227,7 +227,7 @@ namespace MathNet.Numerics
return ans * ax;
}
-
+
///
/// Returns the upper incomplete gamma function
/// Gamma(a,x) = 1/Gamma(a) * int(exp(-t)t^(a-1),t=0..x) for real a > 0, x > 0.
@@ -239,7 +239,7 @@ namespace MathNet.Numerics
{
return GammaUpperRegularized(a, x) * Gamma(a);
}
-
+
///
/// Returns the lower incomplete gamma function
/// gamma(a,x) = int(exp(-t)t^(a-1),t=0..x) for real a > 0, x > 0.
@@ -363,6 +363,195 @@ namespace MathNet.Numerics
return 1d - (Math.Exp(ax) * ans);
}
+ ///
+ /// Returns the inverse P^(-1) of the regularized lower incomplete gamma function
+ /// P(a,x) = 1/Gamma(a) * int(exp(-t)t^(a-1),t=0..x) for real a > 0, x > 0,
+ /// such that P^(-1)(a,P(a,x)) == x.
+ ///
+ public static double GammaLowerRegularizedInv(double a, double y0)
+ {
+ const double epsilon = 0.000000000000001;
+ const double big = 4503599627370496.0;
+ const double threshold = 5*epsilon;
+
+ if (double.IsNaN(a) || double.IsNaN(y0))
+ {
+ return double.NaN;
+ }
+
+ if (a < 0 || a.AlmostEqual(0.0) || y0 < 0 || y0 > 1)
+ {
+ throw new ArgumentOutOfRangeException("a,y0", Properties.Resources.ArgumentNotNegative);
+ }
+
+ if (y0.AlmostEqual(0.0))
+ {
+ return 0d;
+ }
+
+ if (y0.AlmostEqual(1.0))
+ {
+ return Double.PositiveInfinity;
+ }
+
+ y0 = 1 - y0;
+
+ double xUpper = big;
+ double xLower = 0;
+ double yUpper = 1;
+ double yLower = 0;
+
+ // Initial Guess
+ double d = 1/(9*a);
+ double y = 1 - d - (0.98*Constants.Sqrt2*ErfInv((2.0*y0) - 1.0)*Math.Sqrt(d));
+ double x = a*y*y*y;
+ double lgm = GammaLn(a);
+
+ for (int i = 0; i < 10; i++)
+ {
+ if (x < xLower || x > xUpper)
+ {
+ d = 0.0625;
+ break;
+ }
+
+ y = 1 - GammaLowerRegularized(a, x);
+ if (y < yLower || y > yUpper)
+ {
+ d = 0.0625;
+ break;
+ }
+
+ if (y < y0)
+ {
+ xUpper = x;
+ yLower = y;
+ }
+ else
+ {
+ xLower = x;
+ yUpper = y;
+ }
+
+ d = ((a - 1)*Math.Log(x)) - x - lgm;
+ if (d < -709.78271289338399)
+ {
+ d = 0.0625;
+ break;
+ }
+
+ d = -Math.Exp(d);
+ d = (y - y0)/d;
+ if (Math.Abs(d/x) < epsilon)
+ {
+ return x;
+ }
+
+ if ((d > (x/4)) && (y0 < 0.05))
+ {
+ // Naive heuristics for cases near the singularity
+ d = x/10;
+ }
+
+ x -= d;
+ }
+
+ if (xUpper == big)
+ {
+ if (x <= 0)
+ {
+ x = 1;
+ }
+
+ while (xUpper == big)
+ {
+ x = (1 + d)*x;
+ y = 1 - GammaLowerRegularized(a, x);
+ if (y < y0)
+ {
+ xUpper = x;
+ yLower = y;
+ break;
+ }
+
+ d = d + d;
+ }
+ }
+
+ int dir = 0;
+ d = 0.5;
+ for (int i = 0; i < 400; i++)
+ {
+ x = xLower + (d*(xUpper - xLower));
+ y = 1 - GammaLowerRegularized(a, x);
+ lgm = (xUpper - xLower)/(xLower + xUpper);
+ if (Math.Abs(lgm) < threshold)
+ {
+ return x;
+ }
+
+ lgm = (y - y0)/y0;
+ if (Math.Abs(lgm) < threshold)
+ {
+ return x;
+ }
+
+ if (x <= 0d)
+ {
+ return 0d;
+ }
+
+ if (y >= y0)
+ {
+ xLower = x;
+ yUpper = y;
+ if (dir < 0)
+ {
+ dir = 0;
+ d = 0.5;
+ }
+ else
+ {
+ if (dir > 1)
+ {
+ d = (0.5*d) + 0.5;
+ }
+ else
+ {
+ d = (y0 - yLower)/(yUpper - yLower);
+ }
+ }
+
+ dir = dir + 1;
+ }
+ else
+ {
+ xUpper = x;
+ yLower = y;
+ if (dir > 0)
+ {
+ dir = 0;
+ d = 0.5;
+ }
+ else
+ {
+ if (dir < -1)
+ {
+ d = 0.5*d;
+ }
+ else
+ {
+ d = (y0 - yLower)/(yUpper - yLower);
+ }
+ }
+
+ dir = dir - 1;
+ }
+ }
+
+ return x;
+ }
+
///
/// Computes the Digamma function which is mathematically defined as the derivative of the logarithm of the gamma function.
/// This implementation is based on
diff --git a/src/UnitTests/SpecialFunctionsTests/GammaTests.cs b/src/UnitTests/SpecialFunctionsTests/GammaTests.cs
index 714ec108..823f5cb4 100644
--- a/src/UnitTests/SpecialFunctionsTests/GammaTests.cs
+++ b/src/UnitTests/SpecialFunctionsTests/GammaTests.cs
@@ -123,9 +123,36 @@ namespace MathNet.Numerics.UnitTests.SpecialFunctionsTests
[TestCase(1000, 10000, 1.0, 14)]
[TestCase(1e+50, 1e+48, 0.0, 14)]
[TestCase(1e+50, 1e+52, 1.0, 14)]
- public void GammaLowerRegularized(double a, double x, double f, int digits)
+ public void GammaLowerRegularized(double a, double x, double y, int digits)
{
- AssertHelpers.AlmostEqualRelative(f, SpecialFunctions.GammaLowerRegularized(a, x), digits);
+ AssertHelpers.AlmostEqualRelative(y, SpecialFunctions.GammaLowerRegularized(a, x), digits);
+ }
+
+ ///
+ /// Gamma lower regularized inverse.
+ ///
+ [TestCase(double.NaN, Double.NaN, Double.NaN, 14)]
+ [TestCase(0.1, 1.0, 0.97587265627367222115949155252812057714751052498477013, 13)]
+ [TestCase(0.1, 2.0, 0.99432617602018847196075251078067514034772764693462125, 13)]
+ [TestCase(0.1, 8.0, 0.99999507519205198048686442150578226823401842046310854, 10)]
+ [TestCase(1.5, 1.0, 0.42759329552912016600095238564127189392715996802703368, 13)]
+ [TestCase(1.5, 2.0, 0.73853587005088937779717792402407879809718939080920993, 13)]
+ [TestCase(1.5, 8.0, 0.99886601571021467734329986257903021041757398191304284, 13)]
+ [TestCase(2.5, 1.0, 0.15085496391539036377410688601371365034788861473418704, 13)]
+ [TestCase(2.5, 2.0, 0.45058404864721976739416885516693969548484517509263197, 13)]
+ [TestCase(2.5, 8.0, 0.99315592607757956900093935107222761316136944145439676, 13)]
+ [TestCase(5.5, 1.0, 0.0015041182825838038421585211353488839717739161316985392, 13)]
+ [TestCase(5.5, 2.0, 0.030082976121226050615171484772387355162056796585883967, 13)]
+ [TestCase(5.5, 8.0, 0.85886911973294184646060071855669224657735916933487681, 13)]
+ [TestCase(100, 90, 0.1582209891864301681049696996709105316998233457433473, 12)]
+ [TestCase(100, 100, 0.5132987982791486648573142565640291634709251499279450, 12)]
+ [TestCase(100, 110, 0.8417213299399129061982996209829688531933500308658222, 12)]
+ [TestCase(500, 450, 0.0107172380912897415573958770655204965434869949241480, 12)]
+ [TestCase(500, 500, 0.5059471461707603580470479574412058032802735425634263, 12)]
+ [TestCase(500, 550, 0.9853855918737048059548470006900844665580616318702748, 12)]
+ public void GammaLowerRegularizedInv(double a, double x, double y, int digits)
+ {
+ AssertHelpers.AlmostEqualRelative(x, SpecialFunctions.GammaLowerRegularizedInv(a, y), digits);
}
///