Special Functions: inverse of Gamma lower regularized

13 years ago · ad208ae9c7
2 changed files with 224 additions and 8 deletions
--- a/src/Numerics/SpecialFunctions/Gamma.cs
+++ b/src/Numerics/SpecialFunctions/Gamma.cs
@ -70,7 +70,7 @@ namespace MathNet.Numerics
        };

        /// <summary>
-        /// Computes the logarithm of the Gamma function. 
+        /// Computes the logarithm of the Gamma function.
        /// </summary>
        /// <param name="z">The argument of the gamma function.</param>
        /// <returns>The logarithm of the gamma function.</returns>
@ -112,7 +112,7 @@ namespace MathNet.Numerics
        }

        /// <summary>
-        /// Computes the Gamma function. 
+        /// Computes the Gamma function.
        /// </summary>
        /// <param name="z">The argument of the gamma function.</param>
        /// <returns>The logarithm of the gamma function.</returns>
@ -151,7 +151,7 @@ namespace MathNet.Numerics
                return s * Constants.TwoSqrtEOverPi * Math.Pow((z - 0.5 + GammaR) / Math.E, z - 0.5);
            }
        }
-    
+
        /// <summary>
        /// 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 &gt; 0, x &gt; 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;
        }
-    
+
        /// <summary>
        /// Returns the upper incomplete gamma function
        /// Gamma(a,x) = 1/Gamma(a) * int(exp(-t)t^(a-1),t=0..x) for real a &gt; 0, x &gt; 0.
@ -239,7 +239,7 @@ namespace MathNet.Numerics
        {
            return GammaUpperRegularized(a, x) * Gamma(a);
        }
-    
+
        /// <summary>
        /// Returns the lower incomplete gamma function
        /// gamma(a,x) = int(exp(-t)t^(a-1),t=0..x) for real a &gt; 0, x &gt; 0.
@ -363,6 +363,195 @@ namespace MathNet.Numerics
            return 1d - (Math.Exp(ax) * ans);
        }

+        /// <summary>
+        /// 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 &gt; 0, x &gt; 0,
+        /// such that P^(-1)(a,P(a,x)) == x.
+        /// </summary>
+        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;
+        }
+
        /// <summary>
        /// Computes the Digamma function which is mathematically defined as the derivative of the logarithm of the gamma function.
        /// This implementation is based on
--- 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);
+        }
+
+        /// <summary>
+        /// Gamma lower regularized inverse.
+        /// </summary>
+        [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);
        }

        /// <summary>