From e12b1e21626b03b1f1eca3ea35a38ec8e3c280b4 Mon Sep 17 00:00:00 2001 From: Christoph Ruegg Date: Sat, 24 Aug 2013 02:07:17 +0200 Subject: [PATCH] Distributions: adapt Beta --- src/Numerics/Distributions/Beta.cs | 212 ++++++++++-------- .../DistributionTests/Continuous/BetaTests.cs | 3 + 2 files changed, 121 insertions(+), 94 deletions(-) diff --git a/src/Numerics/Distributions/Beta.cs b/src/Numerics/Distributions/Beta.cs index 130b79df..7f27caa1 100644 --- a/src/Numerics/Distributions/Beta.cs +++ b/src/Numerics/Distributions/Beta.cs @@ -63,7 +63,6 @@ namespace MathNet.Numerics.Distributions /// /// The α shape parameter of the Beta distribution. Range: α ≥ 0. /// The β shape parameter of the Beta distribution. Range: β ≥ 0. - /// If any of the Beta parameters are negative. public Beta(double a, double b) { _random = new System.Random(); @@ -76,7 +75,6 @@ namespace MathNet.Numerics.Distributions /// The α shape parameter of the Beta distribution. Range: α ≥ 0. /// The β shape parameter of the Beta distribution. Range: β ≥ 0. /// The random number generator which is used to draw random samples. - /// If any of the Beta parameters are negative. public Beta(double a, double b, System.Random randomSource) { _random = randomSource ?? new System.Random(); @@ -92,17 +90,6 @@ namespace MathNet.Numerics.Distributions return "Beta(α = " + _shapeA + ", β = " + _shapeB + ")"; } - /// - /// Checks whether the parameters of the distribution are valid. - /// - /// The α shape parameter of the Beta distribution. Range: α ≥ 0. - /// The β shape parameter of the Beta distribution. Range: β ≥ 0. - /// true when the parameters are valid, false otherwise. - static bool IsValidParameterSet(double a, double b) - { - return a >= 0.0 && b >= 0.0; - } - /// /// Sets the parameters of the distribution after checking their validity. /// @@ -111,7 +98,7 @@ namespace MathNet.Numerics.Distributions /// When the parameters are out of range. void SetParameters(double a, double b) { - if (Control.CheckDistributionParameters && !IsValidParameterSet(a, b)) + if (a < 0.0 || b < 0.0 || Double.IsNaN(a) || Double.IsNaN(b)) { throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters); } @@ -314,26 +301,99 @@ namespace MathNet.Numerics.Distributions /// /// The location at which to compute the density. /// the density at . + /// public double Density(double x) { + return PDF(_shapeA, _shapeB, x); + } + + /// + /// Computes the log probability density of the distribution (lnPDF) at x, i.e. ln(∂P(X ≤ x)/∂x). + /// + /// The location at which to compute the log density. + /// the log density at . + /// + public double DensityLn(double x) + { + return PDFLn(_shapeA, _shapeB, x); + } + + /// + /// Computes the cumulative distribution (CDF) of the distribution at x, i.e. P(X ≤ x). + /// + /// The location at which to compute the cumulative distribution function. + /// the cumulative distribution at location . + /// + public double CumulativeDistribution(double x) + { + return CDF(_shapeA, _shapeB, x); + } + + /// + /// Generates a sample from the Beta distribution. + /// + /// a sample from the distribution. + public double Sample() + { + return SampleUnchecked(_random, _shapeA, _shapeB); + } + + /// + /// Generates a sequence of samples from the Beta distribution. + /// + /// a sequence of samples from the distribution. + public IEnumerable Samples() + { + while (true) + { + yield return SampleUnchecked(_random, _shapeA, _shapeB); + } + } + + /// + /// Samples Beta distributed random variables by sampling two Gamma variables and normalizing. + /// + /// The random number generator to use. + /// The α shape parameter of the Beta distribution. Range: α ≥ 0. + /// The β shape parameter of the Beta distribution. Range: β ≥ 0. + /// a random number from the Beta distribution. + static double SampleUnchecked(System.Random rnd, double a, double b) + { + var x = Gamma.SampleUnchecked(rnd, a, 1.0); + var y = Gamma.SampleUnchecked(rnd, b, 1.0); + return x / (x + y); + } + + /// + /// Computes the probability density of the distribution (PDF) at x, i.e. ∂P(X ≤ x)/∂x. + /// + /// The α shape parameter of the Beta distribution. Range: α ≥ 0. + /// The β shape parameter of the Beta distribution. Range: β ≥ 0. + /// The location at which to compute the density. + /// the density at . + /// + public static double PDF(double a, double b, double x) + { + if (a < 0.0 || b < 0.0) throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters); + if (x < 0.0 || x > 1.0) return 0.0; - if (Double.IsPositiveInfinity(_shapeA) && Double.IsPositiveInfinity(_shapeB)) + if (Double.IsPositiveInfinity(a) && Double.IsPositiveInfinity(b)) { return x == 0.5 ? Double.PositiveInfinity : 0.0; } - if (Double.IsPositiveInfinity(_shapeA)) + if (Double.IsPositiveInfinity(a)) { return x == 1.0 ? Double.PositiveInfinity : 0.0; } - if (Double.IsPositiveInfinity(_shapeB)) + if (Double.IsPositiveInfinity(b)) { return x == 0.0 ? Double.PositiveInfinity : 0.0; } - if (_shapeA == 0.0 && _shapeB == 0.0) + if (a == 0.0 && b == 0.0) { if (x == 0.0 || x == 1.0) { @@ -343,85 +403,90 @@ namespace MathNet.Numerics.Distributions return 0.0; } - if (_shapeA == 0.0) return x == 0.0 ? Double.PositiveInfinity : 0.0; - if (_shapeB == 0.0) return x == 1.0 ? Double.PositiveInfinity : 0.0; - if (_shapeA == 1.0 && _shapeB == 1.0) return 1.0; + if (a == 0.0) return x == 0.0 ? Double.PositiveInfinity : 0.0; + if (b == 0.0) return x == 1.0 ? Double.PositiveInfinity : 0.0; + if (a == 1.0 && b == 1.0) return 1.0; - var 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); + var bb = SpecialFunctions.Gamma(a + b) / (SpecialFunctions.Gamma(a) * SpecialFunctions.Gamma(b)); + return bb * Math.Pow(x, a - 1.0) * Math.Pow(1.0 - x, b - 1.0); } /// /// Computes the log probability density of the distribution (lnPDF) at x, i.e. ln(∂P(X ≤ x)/∂x). /// - /// The location at which to compute the log density. + /// The α shape parameter of the Beta distribution. Range: α ≥ 0. + /// The β shape parameter of the Beta distribution. Range: β ≥ 0. + /// The location at which to compute the density. /// the log density at . - public double DensityLn(double x) + /// + public static double PDFLn(double a, double b, double x) { + if (a < 0.0 || b < 0.0) throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters); + if (x < 0.0 || x > 1.0) return Double.NegativeInfinity; - if (Double.IsPositiveInfinity(_shapeA) && Double.IsPositiveInfinity(_shapeB)) + if (Double.IsPositiveInfinity(a) && Double.IsPositiveInfinity(b)) { return x == 0.5 ? Double.PositiveInfinity : Double.NegativeInfinity; } - if (Double.IsPositiveInfinity(_shapeA)) + if (Double.IsPositiveInfinity(a)) { return x == 1.0 ? Double.PositiveInfinity : Double.NegativeInfinity; } - if (Double.IsPositiveInfinity(_shapeB)) + if (Double.IsPositiveInfinity(b)) { return x == 0.0 ? Double.PositiveInfinity : Double.NegativeInfinity; } - if (_shapeA == 0.0 && _shapeB == 0.0) + if (a == 0.0 && b == 0.0) { - if (x == 0.0 || x == 1.0) - { - return Double.PositiveInfinity; - } - - return Double.NegativeInfinity; + return x == 0.0 || x == 1.0 ? Double.PositiveInfinity : Double.NegativeInfinity; } - if (_shapeA == 0.0) return x == 0.0 ? Double.PositiveInfinity : Double.NegativeInfinity; - if (_shapeB == 0.0) return x == 1.0 ? Double.PositiveInfinity : Double.NegativeInfinity; - if (_shapeA == 1.0 && _shapeB == 1.0) return 0.0; + if (a == 0.0) return x == 0.0 ? Double.PositiveInfinity : Double.NegativeInfinity; + if (b == 0.0) return x == 1.0 ? Double.PositiveInfinity : Double.NegativeInfinity; + if (a == 1.0 && b == 1.0) return 0.0; - var a = SpecialFunctions.GammaLn(_shapeA + _shapeB) - SpecialFunctions.GammaLn(_shapeA) - SpecialFunctions.GammaLn(_shapeB); - var b = x == 0.0 ? (_shapeA == 1.0 ? 0.0 : Double.NegativeInfinity) : (_shapeA - 1.0)*Math.Log(x); - var c = x == 1.0 ? (_shapeB == 1.0 ? 0.0 : Double.NegativeInfinity) : (_shapeB - 1.0)*Math.Log(1.0 - x); + var aa = SpecialFunctions.GammaLn(a + b) - SpecialFunctions.GammaLn(a) - SpecialFunctions.GammaLn(b); + var bb = x == 0.0 ? (a == 1.0 ? 0.0 : Double.NegativeInfinity) : (a - 1.0)*Math.Log(x); + var cc = x == 1.0 ? (b == 1.0 ? 0.0 : Double.NegativeInfinity) : (b - 1.0)*Math.Log(1.0 - x); - return a + b + c; + return aa + bb + cc; } /// /// Computes the cumulative distribution (CDF) of the distribution at x, i.e. P(X ≤ x). /// /// The location at which to compute the cumulative distribution function. + /// The α shape parameter of the Beta distribution. Range: α ≥ 0. + /// The β shape parameter of the Beta distribution. Range: β ≥ 0.> /// the cumulative distribution at location . - public double CumulativeDistribution(double x) + /// + public static double CDF(double a, double b, double x) { + if (a < 0.0 || b < 0.0) throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters); + if (x < 0.0) return 0.0; if (x >= 1.0) return 1.0; - if (Double.IsPositiveInfinity(_shapeA) && Double.IsPositiveInfinity(_shapeB)) + if (Double.IsPositiveInfinity(a) && Double.IsPositiveInfinity(b)) { return x < 0.5 ? 0.0 : 1.0; } - if (Double.IsPositiveInfinity(_shapeA)) + if (Double.IsPositiveInfinity(a)) { return x < 1.0 ? 0.0 : 1.0; } - if (Double.IsPositiveInfinity(_shapeB)) + if (Double.IsPositiveInfinity(b)) { return x >= 0.0 ? 1.0 : 0.0; } - if (_shapeA == 0.0 && _shapeB == 0.0) + if (a == 0.0 && b == 0.0) { if (x >= 0.0 && x < 1.0) { @@ -431,46 +496,11 @@ namespace MathNet.Numerics.Distributions return 1.0; } - if (_shapeA == 0.0) return 1.0; - if (_shapeB == 0.0) return x >= 1.0 ? 1.0 : 0.0; - if (_shapeA == 1.0 && _shapeB == 1.0) return x; - - return SpecialFunctions.BetaRegularized(_shapeA, _shapeB, x); - } - - /// - /// Samples Beta distributed random variables by sampling two Gamma variables and normalizing. - /// - /// The random number generator to use. - /// The α shape parameter of the Beta distribution. Range: α ≥ 0. - /// The β shape parameter of the Beta distribution. Range: β ≥ 0. - /// a random number from the Beta distribution. - static double SampleUnchecked(System.Random rnd, double a, double b) - { - var x = Gamma.SampleUnchecked(rnd, a, 1.0); - var y = Gamma.SampleUnchecked(rnd, b, 1.0); - return x/(x + y); - } - - /// - /// Generates a sample from the Beta distribution. - /// - /// a sample from the distribution. - public double Sample() - { - return SampleUnchecked(_random, _shapeA, _shapeB); - } + if (a == 0.0) return 1.0; + if (b == 0.0) return x >= 1.0 ? 1.0 : 0.0; + if (a == 1.0 && b == 1.0) return x; - /// - /// Generates a sequence of samples from the Beta distribution. - /// - /// a sequence of samples from the distribution. - public IEnumerable Samples() - { - while (true) - { - yield return SampleUnchecked(_random, _shapeA, _shapeB); - } + return SpecialFunctions.BetaRegularized(a, b, x); } /// @@ -482,10 +512,7 @@ namespace MathNet.Numerics.Distributions /// a sample from the distribution. public static double Sample(System.Random rnd, double a, double b) { - if (Control.CheckDistributionParameters && !IsValidParameterSet(a, b)) - { - throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters); - } + if (a < 0.0 || b < 0.0) throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters); return SampleUnchecked(rnd, a, b); } @@ -499,10 +526,7 @@ namespace MathNet.Numerics.Distributions /// a sequence of samples from the distribution. public static IEnumerable Samples(System.Random rnd, double a, double b) { - if (Control.CheckDistributionParameters && !IsValidParameterSet(a, b)) - { - throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters); - } + if (a < 0.0 || b < 0.0) throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters); while (true) { diff --git a/src/UnitTests/DistributionTests/Continuous/BetaTests.cs b/src/UnitTests/DistributionTests/Continuous/BetaTests.cs index 6381b4b8..7694e583 100644 --- a/src/UnitTests/DistributionTests/Continuous/BetaTests.cs +++ b/src/UnitTests/DistributionTests/Continuous/BetaTests.cs @@ -369,6 +369,7 @@ namespace MathNet.Numerics.UnitTests.DistributionTests.Continuous { var n = new Beta(a, b); AssertHelpers.AlmostEqual(pdf, n.Density(x), 13); + AssertHelpers.AlmostEqual(pdf, Beta.PDF(a, b, x), 13); } /// @@ -414,6 +415,7 @@ namespace MathNet.Numerics.UnitTests.DistributionTests.Continuous { var n = new Beta(a, b); AssertHelpers.AlmostEqual(pdfln, n.DensityLn(x), 14); + AssertHelpers.AlmostEqual(pdfln, Beta.PDFLn(a, b, x), 14); } /// @@ -457,6 +459,7 @@ namespace MathNet.Numerics.UnitTests.DistributionTests.Continuous { var n = new Beta(a, b); AssertHelpers.AlmostEqual(cdf, n.CumulativeDistribution(x), 13); + AssertHelpers.AlmostEqual(cdf, Beta.CDF(a, b, x), 13); } } }