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;
+ }
}
///