From b643907a155a9636693c19f7c111ba372bf53fe7 Mon Sep 17 00:00:00 2001 From: Jurgen Van Gael Date: Mon, 10 Aug 2009 15:30:59 +0800 Subject: [PATCH] Fixed some bugs in the Beta distribution. Signed-off-by: Christoph Ruegg --- .../DistributionTests/Continuous/BetaTests.cs | 14 +- src/Managed/Distributions/Continuous/Beta.cs | 180 +++++++++++++++++- 2 files changed, 184 insertions(+), 10 deletions(-) diff --git a/src/Managed.UnitTests/DistributionTests/Continuous/BetaTests.cs b/src/Managed.UnitTests/DistributionTests/Continuous/BetaTests.cs index 2110c02c..39ce90b3 100644 --- a/src/Managed.UnitTests/DistributionTests/Continuous/BetaTests.cs +++ b/src/Managed.UnitTests/DistributionTests/Continuous/BetaTests.cs @@ -278,8 +278,10 @@ namespace MathNet.Numerics.UnitTests.DistributionTests [Row(1.0, 1.0, 0.5, 1.0)] [Row(1.0, 1.0, 1.0, 1.0)] [Row(9.0, 1.0, 0.0, 0.0)] - [Row(9.0, 1.0, 0.5, 0.035155378090821160189479427593561667617370600927556366)] - [Row(9.0, 1.0, 1.0, 8.9997767912502170085067334639517869100468738374544298)] + [Row(9.0, 1.0, 0.5, 0.03515625)] + [Row(9.0, 1.0, 1.0, 9.0)] + [Row(9.0, 1.0, -1.0, 0.0)] + [Row(9.0, 1.0, 2.0, 0.0)] [Row(5.0, 100.0, 0.0, 0.0)] [Row(5.0, 100.0, 0.5, 1.0881845516040810386311829462908430145307026037926335e-21)] [Row(5.0, 100.0, 1.0, 0.0)] @@ -298,7 +300,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), 14); + AssertHelpers.AlmostEqual(pdf, n.Density(x), 13); } [Test] @@ -315,8 +317,10 @@ namespace MathNet.Numerics.UnitTests.DistributionTests [Row(1.0, 1.0, 0.5, 0.0)] [Row(1.0, 1.0, 1.0, 0.0)] [Row(9.0, 1.0, 0.0, Double.NegativeInfinity)] - [Row(9.0, 1.0, 0.5, -3.3479528671433430925473664978203611353090199592365404)] - [Row(9.0, 1.0, 1.0, 2.197199776056471990627201569939718235188380979206923)] + [Row(9.0, 1.0, 0.5, -3.3479528671433430925473664978203611353090199592365458)] + [Row(9.0, 1.0, 1.0, 2.1972245773362193827904904738450514092949811156454996)] + [Row(9.0, 1.0, -1.0, Double.NegativeInfinity)] + [Row(9.0, 1.0, 2.0, Double.NegativeInfinity)] [Row(5.0, 100.0, 0.0, Double.NegativeInfinity)] [Row(5.0, 100.0, 0.5, -51.447830024537682154565870837960406410586196074573801)] [Row(5.0, 100.0, 1.0, Double.NegativeInfinity)] diff --git a/src/Managed/Distributions/Continuous/Beta.cs b/src/Managed/Distributions/Continuous/Beta.cs index 7433793b..31773bfe 100644 --- a/src/Managed/Distributions/Continuous/Beta.cs +++ b/src/Managed/Distributions/Continuous/Beta.cs @@ -236,12 +236,24 @@ namespace MathNet.Numerics.Distributions } else if (Double.IsPositiveInfinity(_shapeA)) { - return 2.0; + return -2.0; } else if (Double.IsPositiveInfinity(_shapeB)) { return 2.0; } + else if (_shapeA == 0.0 && _shapeB == 0.0) + { + return 0.0; + } + else if (_shapeA == 0.0) + { + return 2.0; + } + else if (_shapeB == 0.0) + { + return -2.0; + } else { return 2.0 * (_shapeB - _shapeA) * Math.Sqrt(_shapeA + _shapeB + 1.0) @@ -326,8 +338,86 @@ namespace MathNet.Numerics.Distributions /// the density at . public double Density(double x) { - double b = SpecialFunctions.Gamma(_shapeA + _shapeB) / (SpecialFunctions.Gamma(_shapeA) * SpecialFunctions.Gamma(_shapeB)); - return b * Math.Pow(x, _shapeA - 1.0) * Math.Pow(1.0 - x, _shapeB - 1.0); + if (x < 0.0 || x > 1.0) + { + return 0.0; + } + + if (Double.IsPositiveInfinity(_shapeA) && Double.IsPositiveInfinity(_shapeB)) + { + if (x == 0.5) + { + return Double.PositiveInfinity; + } + else + { + return 0.0; + } + } + else if (Double.IsPositiveInfinity(_shapeA)) + { + if (x == 1.0) + { + return Double.PositiveInfinity; + } + else + { + return 0.0; + } + } + else if (Double.IsPositiveInfinity(_shapeB)) + { + if (x == 0.0) + { + return Double.PositiveInfinity; + } + else + { + return 0.0; + } + } + else if (_shapeA == 0.0 && _shapeB == 0.0) + { + if (x == 0.0 || x == 1.0) + { + return Double.PositiveInfinity; + } + else + { + return 0.0; + } + } + else if (_shapeA == 0.0) + { + if (x == 0.0) + { + return Double.PositiveInfinity; + } + else + { + return 0.0; + } + } + else if (_shapeB == 0.0) + { + if (x == 1.0) + { + return Double.PositiveInfinity; + } + else + { + return 0.0; + } + } + else if (_shapeA == 1.0 && _shapeB == 1.0) + { + return 1.0; + } + else + { + double b = SpecialFunctions.Gamma(_shapeA + _shapeB) / (SpecialFunctions.Gamma(_shapeA) * SpecialFunctions.Gamma(_shapeB)); + return b * Math.Pow(x, _shapeA - 1.0) * Math.Pow(1.0 - x, _shapeB - 1.0); + } } /// @@ -337,8 +427,88 @@ namespace MathNet.Numerics.Distributions /// the log density at . public double DensityLn(double x) { - double b = SpecialFunctions.GammaLn(_shapeA + _shapeB) - SpecialFunctions.GammaLn(_shapeA) - SpecialFunctions.GammaLn(_shapeB); - return b + (_shapeA - 1.0)*Math.Log(x) + (_shapeB - 1.0)*Math.Log(1.0 - x); + if (x < 0.0 || x > 1.0) + { + return Double.NegativeInfinity; + } + + if (Double.IsPositiveInfinity(_shapeA) && Double.IsPositiveInfinity(_shapeB)) + { + if (x == 0.5) + { + return Double.PositiveInfinity; + } + else + { + return Double.NegativeInfinity; + } + } + else if (Double.IsPositiveInfinity(_shapeA)) + { + if (x == 1.0) + { + return Double.PositiveInfinity; + } + else + { + return Double.NegativeInfinity; + } + } + else if (Double.IsPositiveInfinity(_shapeB)) + { + if (x == 0.0) + { + return Double.PositiveInfinity; + } + else + { + return Double.NegativeInfinity; + } + } + else if (_shapeA == 0.0 && _shapeB == 0.0) + { + if (x == 0.0 || x == 1.0) + { + return Double.PositiveInfinity; + } + else + { + return Double.NegativeInfinity; + } + } + else if (_shapeA == 0.0) + { + if (x == 0.0) + { + return Double.PositiveInfinity; + } + else + { + return Double.NegativeInfinity; + } + } + else if (_shapeB == 0.0) + { + if (x == 1.0) + { + return Double.PositiveInfinity; + } + else + { + return Double.NegativeInfinity; + } + } + else if (_shapeA == 1.0 && _shapeB == 1.0) + { + return 0.0; + } + else + { + double a = SpecialFunctions.GammaLn(_shapeA + _shapeB) - SpecialFunctions.GammaLn(_shapeA) - SpecialFunctions.GammaLn(_shapeB); + double b = (x == 0.0 ? (_shapeA == 1.0 ? 0.0 : Double.NegativeInfinity) : (_shapeA - 1.0) * Math.Log(x)); + double c = (x == 1.0 ? (_shapeB == 1.0 ? 0.0 : Double.NegativeInfinity) : (_shapeB - 1.0) * Math.Log(1.0 - x)); + return a + b + c; + } } ///