Added unit tests for Gamma & Beta.

Added DiGamma function.

Signed-off-by: Christoph Ruegg <git@cdrnet.ch>

src/Managed.UnitTests/DistributionTests/Continuous/BetaTests.cs

@@ -170,7 +170,7 @@ namespace MathNet.Numerics.UnitTests.DistributionTests
         public void ValidateSkewness(double a, double b, double skewness)
         {
             var n = new Beta(a, b);
-            AssertEx.AreEqual<double>(skewness, n.Skewness);
+            AssertHelpers.AlmostEqual(skewness, n.Skewness, 15);
         }
 
         [Test]
@@ -298,7 +298,7 @@ namespace MathNet.Numerics.UnitTests.DistributionTests
         public void ValidateDensity(double a, double b, double x, double pdf)
         {
             var n = new Beta(a, b);
-            AssertHelpers.AlmostEqual(pdf, n.Density(x), 15);
+            AssertHelpers.AlmostEqual(pdf, n.Density(x), 14);
         }
 
         [Test]
@@ -335,7 +335,7 @@ namespace MathNet.Numerics.UnitTests.DistributionTests
         public void ValidateDensityLn(double a, double b, double x, double pdfln)
         {
             var n = new Beta(a, b);
-            AssertHelpers.AlmostEqual(pdfln, n.DensityLn(x), 15);
+            AssertHelpers.AlmostEqual(pdfln, n.DensityLn(x), 14);
         }
 
         [Test, Ignore("Depending on Special Functions")]

src/Managed.UnitTests/DistributionTests/Continuous/GammaTests.cs

@@ -181,7 +181,7 @@ namespace MathNet.Numerics.UnitTests.DistributionTests
         }
 
         [Test]
-        [Row(0.0, 0.0, 0.0)]
+        [Row(0.0, 0.0, Double.NaN)]
         [Row(1.0, 0.1, 100.0)]
         [Row(1.0, 1.0, 1.0)]
         [Row(10.0, 10.0, 0.1)]
@@ -190,11 +190,11 @@ namespace MathNet.Numerics.UnitTests.DistributionTests
         public void ValidateVariance(double shape, double invScale, double var)
         {
             var n = new Gamma(shape, invScale);
-            AssertEx.AreEqual<double>(var, n.Variance);
+            AssertHelpers.AlmostEqual(var, n.Variance, 15);
         }
 
         [Test]
-        [Row(0.0, 0.0, 0.0)]
+        [Row(0.0, 0.0, Double.NaN)]
         [Row(1.0, 0.1, 10.0)]
         [Row(1.0, 1.0, 1.0)]
         [Row(10.0, 10.0, 0.31622776601683794197697302588502426416723164097476643)]
@@ -206,8 +206,8 @@ namespace MathNet.Numerics.UnitTests.DistributionTests
             AssertHelpers.AlmostEqual(sdev, n.StdDev, 15);
         }
 
-        [Test, Ignore("Depending on Special Functions")]
-        [Row(0.0, 0.0, Double.PositiveInfinity)]
+        [Test]
+        [Row(0.0, 0.0, Double.NaN)]
         [Row(1.0, 0.1, 3.3025850929940456285068402234265387271634735938763824)]
         [Row(1.0, 1.0, 1.0)]
         [Row(10.0, 10.0, 0.23346908548693395836262094490967812177376750477943892)]
@@ -216,11 +216,11 @@ namespace MathNet.Numerics.UnitTests.DistributionTests
         public void ValidateEntropy(double shape, double invScale, double entropy)
         {
             var n = new Gamma(shape, invScale);
-            AssertHelpers.AlmostEqual(entropy, n.Entropy, 15);
+            AssertHelpers.AlmostEqual(entropy, n.Entropy, 13);
         }
 
         [Test]
-        [Row(0.0, 0.0, Double.PositiveInfinity)]
+        [Row(0.0, 0.0, Double.NaN)]
         [Row(1.0, 0.1, 2.0)]
         [Row(1.0, 1.0, 2.0)]
         [Row(10.0, 10.0, 0.63245553203367586639977870888654370674391102786504337)]
@@ -233,7 +233,7 @@ namespace MathNet.Numerics.UnitTests.DistributionTests
         }
 
         [Test]
-        [Row(0.0, 0.0, Double.PositiveInfinity)]
+        [Row(0.0, 0.0, Double.NaN)]
         [Row(1.0, 0.1, 0.0)]
         [Row(1.0, 1.0, 0.0)]
         [Row(10.0, 10.0, 0.9)]
@@ -253,7 +253,7 @@ namespace MathNet.Numerics.UnitTests.DistributionTests
         [Row(10.0, 10.0)]
         [Row(10.0, 1.0)]
         [Row(10.0, Double.PositiveInfinity)]
-        public void ValidateMedian(double shape, double invScale, double mode)
+        public void ValidateMedian(double shape, double invScale)
         {
             var n = new Gamma(shape, invScale);
             var median = n.Median;

src/Managed.UnitTests/SpecialFunctionsTest/SpecialFunctionsTests.cs

@@ -77,5 +77,28 @@ namespace MathNet.Numerics.UnitTests.SpecialFunctionTests
         {
             AssertHelpers.AlmostEqual(f, SpecialFunctions.Gamma(x), 13);
         }
+
+        [Test]
+        [Row(Double.NaN, Double.NaN)]
+        [Row(-1.5, 0.70315664064524318722569033366791109947350706200623256)]
+        [Row(-0.5, 0.036489973978576520559023667001244432806840395339565891)]
+        [Row(0.1, -10.423754940411076232100295314502760886768558023951363)]
+        [Row(1.0, -0.57721566490153286060651209008240243104215933593992359)]
+        [Row(1.5, 0.036489973978576520559023667001244432806840395339565888)]
+        [Row(Constants.Pi / 2, 0.10067337642740238636795561404029690452798358068944001)]
+        [Row(2.0, 0.42278433509846713939348790991759756895784066406007641)]
+        [Row(2.5, 0.70315664064524318722569033366791109947350706200623255)]
+        [Row(3.0, 0.92278433509846713939348790991759756895784066406007641)]
+        [Row(Constants.Pi, 0.97721330794200673329206948640618234364083460999432603)]
+        [Row(3.5, 1.1031566406452431872256903336679110994735070620062326)]
+        [Row(4.0, 1.2561176684318004727268212432509309022911739973934097)]
+        [Row(4.5, 1.3888709263595289015114046193821968137592213477205183)]
+        [Row(5.0, 1.5061176684318004727268212432509309022911739973934097)]
+        [Row(5.5, 1.6110931485817511237336268416044190359814435699427405)]
+        [Row(10.1, 2.2622143570941481235561593642219403924532310597356171)]
+        public void DiGamma(double x, double f)
+        {
+            AssertHelpers.AlmostEqual(f, SpecialFunctions.DiGamma(x), 13);
+        }
     }
 }

src/Managed/Distributions/Continuous/Beta.cs

@@ -230,8 +230,23 @@ namespace MathNet.Numerics.Distributions
         {
             get
             {
-                return 2.0 * (_shapeB - _shapeA) * Math.Sqrt(_shapeA + _shapeB + 1.0) 
-                    / ((_shapeA + _shapeB + 2.0) * Math.Sqrt(_shapeA * _shapeB));
+                if (Double.IsPositiveInfinity(_shapeA) && Double.IsPositiveInfinity(_shapeB))
+                {
+                    return 0.0;
+                }
+                else if (Double.IsPositiveInfinity(_shapeA))
+                {
+                    return 2.0;
+                }
+                else if (Double.IsPositiveInfinity(_shapeB))
+                {
+                    return 2.0;
+                }
+                else
+                {
+                    return 2.0 * (_shapeB - _shapeA) * Math.Sqrt(_shapeA + _shapeB + 1.0)
+                       / ((_shapeA + _shapeB + 2.0) * Math.Sqrt(_shapeA * _shapeB));
+                }
             }
         }
         #endregion

src/Managed/Distributions/Continuous/Gamma.cs

@@ -39,7 +39,8 @@ namespace MathNet.Numerics.Distributions
     /// <remarks>
     /// <para>The Gamma distribution is parametrized by a shape and inverse scale parameter. When we want
     /// to specify a Gamma distribution which is a point distribution we set the shape parameter to be the
-    /// location of the point distribution and the inverse scale as positive infinity.</para>
+    /// location of the point distribution and the inverse scale as positive infinity. The distribution
+    /// with shape and inverse scale both zero is undefined.</para>
     /// <para> Random number generation for the Gamma distribution is based on the algorithm in:
     /// "A Simple Method for Generating Gamma Variables" - Marsaglia &amp; Tsang
     /// ACM Transactions on Mathematical Software, Vol. 26, No. 3, September 2000, Pages 363–372.</para>
@@ -253,6 +254,10 @@ namespace MathNet.Numerics.Distributions
                 {
                     return 0.0;
                 }
+                else if (_invScale == 0.0 && _shape == 0.0)
+                {
+                    return Double.NaN;
+                }
                 else
                 {
                     return _shape / (_invScale * _invScale);
@@ -271,6 +276,10 @@ namespace MathNet.Numerics.Distributions
                 {
                     return 0.0;
                 }
+                else if (_invScale == 0.0 && _shape == 0.0)
+                {
+                    return Double.NaN;
+                }
                 else
                 {
                     return Math.Sqrt(_shape / (_invScale * _invScale));
@@ -289,6 +298,10 @@ namespace MathNet.Numerics.Distributions
                 {
                     return 0.0;
                 }
+                else if (_invScale == 0.0 && _shape == 0.0)
+                {
+                    return Double.NaN;
+                }
                 else
                 {
                     return _shape - Math.Log(_invScale) + SpecialFunctions.GammaLn(_shape) + (1.0 - _shape) * SpecialFunctions.DiGamma(_shape);
@@ -301,7 +314,21 @@ namespace MathNet.Numerics.Distributions
         /// </summary>
         public double Skewness
         {
-            get { return 2.0 / Math.Sqrt(_shape); }
+            get
+            {
+                if (Double.IsPositiveInfinity(_invScale))
+                {
+                    return 0.0;
+                }
+                else if (_invScale == 0.0 && _shape == 0.0)
+                {
+                    return Double.NaN;
+                }
+                else
+                {
+                    return 2.0 / Math.Sqrt(_shape);
+                }
+            }
         }
         #endregion
 
@@ -318,6 +345,10 @@ namespace MathNet.Numerics.Distributions
                 {
                     return _shape;
                 }
+                else if (_invScale == 0.0 && _shape == 0.0)
+                {
+                    return Double.NaN;
+                }
                 else
                 {
                     return (_shape - 1.0) / _invScale;

src/Managed/SpecialFunctions.cs

@@ -157,6 +157,71 @@ namespace MathNet.Numerics
             }
         }
 
+        /// <summary>
+        /// Computes the digamma function which is mathematically defined as the derivative of the logarithm of the gamma function.
+        /// This implementation is based on
+        ///     Jose Bernardo
+        ///     Algorithm AS 103:
+        ///     Psi ( Digamma ) Function,
+        ///     Applied Statistics,
+        ///     Volume 25, Number 3, 1976, pages 315-317.
+        /// Using the modifications as in Tom Minka's lightspeed toolbox.
+        /// </summary>
+        /// <param name="x">The argument of the digamma function.</param>
+        /// <returns>The value of the DiGamma function at <paramref name="x"/>.</returns>
+        static public double DiGamma(double x)
+        {
+            const double c = 12.0,
+                d1 = -0.57721566490153286,
+                d2 = 1.6449340668482264365,
+                s = 1e-6,
+                s3 = 1.0 / 12.0,
+                s4 = 1.0 / 120.0,
+                s5 = 1.0 / 252.0,
+                s6 = 1.0 / 240.0,
+                s7 = 1.0 / 132.0;
+
+            if (System.Double.IsNegativeInfinity(x) || System.Double.IsNaN(x))
+            {
+                return System.Double.NaN;
+            }
+
+            // Handle special cases.
+            if (x <= 0 && System.Math.Floor(x) == x)
+            {
+                return System.Double.NegativeInfinity;
+            }
+
+            // Use inversion formula for negative numbers.
+            if (x < 0)
+            {
+                return DiGamma(1.0 - x) + System.Math.PI / System.Math.Tan(-System.Math.PI * x);
+            }
+
+            if (x <= s)
+            {
+                return d1 - 1 / x + d2 * x;
+            }
+
+            double result = 0;
+            while (x < c)
+            {
+                result -= 1 / x;
+                x++;
+            }
+
+            if (x >= c)
+            {
+                double r = 1 / x;
+                result += System.Math.Log(x) - 0.5 * r;
+                r *= r;
+
+                result -= r * (s3 - r * (s4 - r * (s5 - r * (s6 - r * s7))));
+            }
+
+            return result;
+        }
+
         public static double IncompleteGamma(double x, double z, bool reg)
         {
             throw new NotImplementedException();
@@ -170,10 +235,5 @@ namespace MathNet.Numerics
         {
             throw new NotImplementedException();
         }
-
-        public static double DiGamma(double x)
-        {
-            throw new NotImplementedException();
-        }
     }
 }