From e9adc7e563ea5be266a1bf7c612d198d08a42c0f Mon Sep 17 00:00:00 2001 From: Christoph Ruegg Date: Sat, 24 Aug 2013 12:52:01 +0200 Subject: [PATCH] Distributions: adapt Gamma --- src/Numerics/Distributions/Gamma.cs | 166 +++++++++--------- .../Continuous/GammaTests.cs | 3 + 2 files changed, 87 insertions(+), 82 deletions(-) diff --git a/src/Numerics/Distributions/Gamma.cs b/src/Numerics/Distributions/Gamma.cs index d1ab96f1..fd2fef04 100644 --- a/src/Numerics/Distributions/Gamma.cs +++ b/src/Numerics/Distributions/Gamma.cs @@ -115,17 +115,6 @@ namespace MathNet.Numerics.Distributions return "Gamma(α = " + _shape + ", β = " + _rate + ")"; } - /// - /// Checks whether the parameters of the distribution are valid. - /// - /// The shape (k, α) of the Gamma distribution. Range: α ≥ 0. - /// The rate or inverse scale (β) of the Gamma distribution. Range: β ≥ 0. - /// true when the parameters are valid, false otherwise. - static bool IsValidParameterSet(double shape, double rate) - { - return shape >= 0.0 && rate >= 0.0; - } - /// /// Sets the parameters of the distribution after checking their validity. /// @@ -134,7 +123,7 @@ namespace MathNet.Numerics.Distributions /// When the parameters are out of range. void SetParameters(double shape, double rate) { - if (Control.CheckDistributionParameters && !IsValidParameterSet(shape, rate)) + if (shape < 0.0 || rate < 0.0 || Double.IsNaN(shape) || Double.IsNaN(rate)) { throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters); } @@ -342,24 +331,10 @@ namespace MathNet.Numerics.Distributions /// /// The location at which to compute the density. /// the density at . + /// public double Density(double x) { - if (Double.IsPositiveInfinity(_rate)) - { - return x == _shape ? Double.PositiveInfinity : 0.0; - } - - if (_shape == 0.0 && _rate == 0.0) - { - return 0.0; - } - - if (_shape == 1.0) - { - return _rate*Math.Exp(-_rate*x); - } - - return Math.Pow(_rate, _shape)*Math.Pow(x, _shape - 1.0)*Math.Exp(-_rate*x)/SpecialFunctions.Gamma(_shape); + return PDF(_shape, _rate, x); } /// @@ -367,24 +342,10 @@ namespace MathNet.Numerics.Distributions /// /// The location at which to compute the log density. /// the log density at . + /// public double DensityLn(double x) { - if (Double.IsPositiveInfinity(_rate)) - { - return x == _shape ? Double.PositiveInfinity : Double.NegativeInfinity; - } - - if (_shape == 0.0 && _rate == 0.0) - { - return Double.NegativeInfinity; - } - - if (_shape == 1.0) - { - return Math.Log(_rate) - (_rate*x); - } - - return (_shape*Math.Log(_rate)) + ((_shape - 1.0)*Math.Log(x)) - (_rate*x) - SpecialFunctions.GammaLn(_shape); + return PDFLn(_shape, _rate, x); } /// @@ -392,19 +353,31 @@ namespace MathNet.Numerics.Distributions /// /// The location at which to compute the cumulative distribution function. /// the cumulative distribution at location . + /// public double CumulativeDistribution(double x) { - if (Double.IsPositiveInfinity(_rate)) - { - return x >= _shape ? 1.0 : 0.0; - } + return CDF(_shape, _rate, x); + } + + /// + /// Generates a sample from the Gamma distribution. + /// + /// a sample from the distribution. + public double Sample() + { + return SampleUnchecked(_random, _shape, _rate); + } - if (_shape == 0.0 && _rate == 0.0) + /// + /// Generates a sequence of samples from the Gamma distribution. + /// + /// a sequence of samples from the distribution. + public IEnumerable Samples() + { + while (true) { - return 0.0; + yield return SampleUnchecked(_random, _shape, _rate); } - - return SpecialFunctions.GammaLowerRegularized(_shape, x*_rate); } /// @@ -431,55 +404,90 @@ namespace MathNet.Numerics.Distributions if (shape < 1.0) { a = shape + 1.0; - alphafix = Math.Pow(rnd.NextDouble(), 1.0/shape); + alphafix = Math.Pow(rnd.NextDouble(), 1.0 / shape); } - var d = a - (1.0/3.0); - var c = 1.0/Math.Sqrt(9.0*d); + var d = a - (1.0 / 3.0); + var c = 1.0 / Math.Sqrt(9.0 * d); while (true) { var x = Normal.Sample(rnd, 0.0, 1.0); - var v = 1.0 + (c*x); + var v = 1.0 + (c * x); while (v <= 0.0) { x = Normal.Sample(rnd, 0.0, 1.0); - v = 1.0 + (c*x); + v = 1.0 + (c * x); } - v = v*v*v; + v = v * v * v; var u = rnd.NextDouble(); - x = x*x; - if (u < 1.0 - (0.0331*x*x)) + x = x * x; + if (u < 1.0 - (0.0331 * x * x)) { - return alphafix*d*v/rate; + return alphafix * d * v / rate; } - if (Math.Log(u) < (0.5*x) + (d*(1.0 - v + Math.Log(v)))) + if (Math.Log(u) < (0.5 * x) + (d * (1.0 - v + Math.Log(v)))) { - return alphafix*d*v/rate; + return alphafix * d * v / rate; } } } /// - /// Generates a sample from the Gamma distribution. + /// Computes the probability density of the distribution (PDF) at x, i.e. ∂P(X ≤ x)/∂x. /// - /// a sample from the distribution. - public double Sample() + /// The shape (k, α) of the Gamma distribution. Range: α ≥ 0. + /// The rate or inverse scale (β) of the Gamma distribution. Range: β ≥ 0. + /// The location at which to compute the density. + /// the density at . + /// + public static double PDF(double shape, double rate, double x) { - return SampleUnchecked(_random, _shape, _rate); + if (shape < 0.0 || rate < 0.0) throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters); + + if (Double.IsPositiveInfinity(rate)) return x == shape ? Double.PositiveInfinity : 0.0; + if (shape == 0.0 && rate == 0.0) return 0.0; + if (shape == 1.0) return rate*Math.Exp(-rate*x); + + return Math.Pow(rate, shape)*Math.Pow(x, shape - 1.0)*Math.Exp(-rate*x)/SpecialFunctions.Gamma(shape); } /// - /// Generates a sequence of samples from the Gamma distribution. + /// Computes the log probability density of the distribution (lnPDF) at x, i.e. ln(∂P(X ≤ x)/∂x). /// - /// a sequence of samples from the distribution. - public IEnumerable Samples() + /// The shape (k, α) of the Gamma distribution. Range: α ≥ 0. + /// The rate or inverse scale (β) of the Gamma distribution. Range: β ≥ 0. + /// The location at which to compute the density. + /// the log density at . + /// + public static double PDFLn(double shape, double rate, double x) { - while (true) - { - yield return SampleUnchecked(_random, _shape, _rate); - } + if (shape < 0.0 || rate < 0.0) throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters); + + if (Double.IsPositiveInfinity(rate)) return x == shape ? Double.PositiveInfinity : Double.NegativeInfinity; + if (shape == 0.0 && rate == 0.0) return Double.NegativeInfinity; + if (shape == 1.0) return Math.Log(rate) - (rate*x); + + return (shape*Math.Log(rate)) + ((shape - 1.0)*Math.Log(x)) - (rate*x) - SpecialFunctions.GammaLn(shape); + } + + /// + /// 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 (k, α) of the Gamma distribution. Range: α ≥ 0. + /// The rate or inverse scale (β) of the Gamma distribution. Range: β ≥ 0. + /// the cumulative distribution at location . + /// + public static double CDF(double shape, double rate, double x) + { + if (shape < 0.0 || rate < 0.0) throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters); + + if (Double.IsPositiveInfinity(rate)) return x >= shape ? 1.0 : 0.0; + if (shape == 0.0 && rate == 0.0) return 0.0; + + return SpecialFunctions.GammaLowerRegularized(shape, x*rate); } /// @@ -491,10 +499,7 @@ namespace MathNet.Numerics.Distributions /// a sample from the distribution. public static double Sample(System.Random rnd, double shape, double rate) { - if (Control.CheckDistributionParameters && !IsValidParameterSet(shape, rate)) - { - throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters); - } + if (shape < 0.0 || rate < 0.0) throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters); return SampleUnchecked(rnd, shape, rate); } @@ -508,10 +513,7 @@ namespace MathNet.Numerics.Distributions /// a sequence of samples from the distribution. public static IEnumerable Samples(System.Random rnd, double shape, double rate) { - if (Control.CheckDistributionParameters && !IsValidParameterSet(shape, rate)) - { - throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters); - } + if (shape < 0.0 || rate < 0.0) throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters); while (true) { diff --git a/src/UnitTests/DistributionTests/Continuous/GammaTests.cs b/src/UnitTests/DistributionTests/Continuous/GammaTests.cs index a7e6b2bf..6442a3d1 100644 --- a/src/UnitTests/DistributionTests/Continuous/GammaTests.cs +++ b/src/UnitTests/DistributionTests/Continuous/GammaTests.cs @@ -378,6 +378,7 @@ namespace MathNet.Numerics.UnitTests.DistributionTests.Continuous { var n = new Gamma(shape, invScale); AssertHelpers.AlmostEqual(pdf, n.Density(x), 14); + AssertHelpers.AlmostEqual(pdf, Gamma.PDF(shape, invScale, x), 14); } /// @@ -409,6 +410,7 @@ namespace MathNet.Numerics.UnitTests.DistributionTests.Continuous { var n = new Gamma(shape, invScale); AssertHelpers.AlmostEqual(pdfln, n.DensityLn(x), 14); + AssertHelpers.AlmostEqual(pdfln, Gamma.PDFLn(shape, invScale, x), 14); } /// @@ -498,6 +500,7 @@ namespace MathNet.Numerics.UnitTests.DistributionTests.Continuous { var n = new Gamma(shape, invScale); AssertHelpers.AlmostEqual(cdf, n.CumulativeDistribution(x), 14); + AssertHelpers.AlmostEqual(cdf, Gamma.CDF(shape, invScale, x), 14); } } }