From 58a2e8693068e934538199dc7626176de00f788d Mon Sep 17 00:00:00 2001 From: jvangael Date: Mon, 10 Aug 2009 14:18:53 +0800 Subject: [PATCH] Added unit tests for Gamma & Beta. Added DiGamma function. Signed-off-by: Christoph Ruegg --- .../DistributionTests/Continuous/BetaTests.cs | 6 +- .../Continuous/GammaTests.cs | 18 ++--- .../SpecialFunctionsTests.cs | 23 ++++++ src/Managed/Distributions/Continuous/Beta.cs | 19 ++++- src/Managed/Distributions/Continuous/Gamma.cs | 35 +++++++++- src/Managed/SpecialFunctions.cs | 70 +++++++++++++++++-- 6 files changed, 150 insertions(+), 21 deletions(-) diff --git a/src/Managed.UnitTests/DistributionTests/Continuous/BetaTests.cs b/src/Managed.UnitTests/DistributionTests/Continuous/BetaTests.cs index abab4984..2110c02c 100644 --- a/src/Managed.UnitTests/DistributionTests/Continuous/BetaTests.cs +++ b/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(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")] diff --git a/src/Managed.UnitTests/DistributionTests/Continuous/GammaTests.cs b/src/Managed.UnitTests/DistributionTests/Continuous/GammaTests.cs index fec6c4a2..5911ce32 100644 --- a/src/Managed.UnitTests/DistributionTests/Continuous/GammaTests.cs +++ b/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(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; diff --git a/src/Managed.UnitTests/SpecialFunctionsTest/SpecialFunctionsTests.cs b/src/Managed.UnitTests/SpecialFunctionsTest/SpecialFunctionsTests.cs index 2f5cf2e1..3e7dd605 100644 --- a/src/Managed.UnitTests/SpecialFunctionsTest/SpecialFunctionsTests.cs +++ b/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); + } } } diff --git a/src/Managed/Distributions/Continuous/Beta.cs b/src/Managed/Distributions/Continuous/Beta.cs index 0d62e2a1..7433793b 100644 --- a/src/Managed/Distributions/Continuous/Beta.cs +++ b/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 diff --git a/src/Managed/Distributions/Continuous/Gamma.cs b/src/Managed/Distributions/Continuous/Gamma.cs index ba103dde..ea4cdb07 100644 --- a/src/Managed/Distributions/Continuous/Gamma.cs +++ b/src/Managed/Distributions/Continuous/Gamma.cs @@ -39,7 +39,8 @@ namespace MathNet.Numerics.Distributions /// /// 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. + /// location of the point distribution and the inverse scale as positive infinity. The distribution + /// with shape and inverse scale both zero is undefined. /// Random number generation for the Gamma distribution is based on the algorithm in: /// "A Simple Method for Generating Gamma Variables" - Marsaglia & Tsang /// ACM Transactions on Mathematical Software, Vol. 26, No. 3, September 2000, Pages 363–372. @@ -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 /// 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; diff --git a/src/Managed/SpecialFunctions.cs b/src/Managed/SpecialFunctions.cs index 858d1947..4761e779 100644 --- a/src/Managed/SpecialFunctions.cs +++ b/src/Managed/SpecialFunctions.cs @@ -157,6 +157,71 @@ namespace MathNet.Numerics } } + /// + /// 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. + /// + /// The argument of the digamma function. + /// The value of the DiGamma function at . + 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(); - } } }