diff --git a/src/Numerics/Distributions/Erlang.cs b/src/Numerics/Distributions/Erlang.cs index eed883b9..d7e30d67 100644 --- a/src/Numerics/Distributions/Erlang.cs +++ b/src/Numerics/Distributions/Erlang.cs @@ -106,17 +106,6 @@ namespace MathNet.Numerics.Distributions return "Erlang(k = " + _shape + ", λ = " + _rate + ")"; } - /// - /// Checks whether the parameters of the distribution are valid. - /// - /// The shape (k) of the Erlang distribution. Range: k ≥ 0. - /// The rate or inverse scale (λ) of the Erlang 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. /// @@ -125,7 +114,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); } @@ -338,24 +327,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); } /// @@ -363,24 +338,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); } /// @@ -388,19 +349,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 Erlang 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 Erlang 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); } /// @@ -427,55 +400,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 Erlang 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 Erlang distribution. Range: k ≥ 0. + /// The rate or inverse scale (λ) of the Erlang 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 Erlang 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 Erlang distribution. Range: k ≥ 0. + /// The rate or inverse scale (λ) of the Erlang 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 Erlang distribution. Range: k ≥ 0. + /// The rate or inverse scale (λ) of the Erlang 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); } /// @@ -487,10 +495,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); } @@ -504,10 +509,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/ErlangTests.cs b/src/UnitTests/DistributionTests/Continuous/ErlangTests.cs index 786b4e6a..60b52b21 100644 --- a/src/UnitTests/DistributionTests/Continuous/ErlangTests.cs +++ b/src/UnitTests/DistributionTests/Continuous/ErlangTests.cs @@ -373,6 +373,7 @@ namespace MathNet.Numerics.UnitTests.DistributionTests.Continuous { var n = new Erlang(shape, invScale); AssertHelpers.AlmostEqual(pdf, n.Density(x), 14); + AssertHelpers.AlmostEqual(pdf, Erlang.PDF(shape, invScale, x), 14); } /// @@ -404,6 +405,7 @@ namespace MathNet.Numerics.UnitTests.DistributionTests.Continuous { var n = new Erlang(shape, invScale); AssertHelpers.AlmostEqual(pdfln, n.DensityLn(x), 14); + AssertHelpers.AlmostEqual(pdfln, Erlang.PDFLn(shape, invScale, x), 14); } /// @@ -456,6 +458,7 @@ namespace MathNet.Numerics.UnitTests.DistributionTests.Continuous { var n = new Erlang(shape, invScale); AssertHelpers.AlmostEqual(cdf, n.CumulativeDistribution(x), 14); + AssertHelpers.AlmostEqual(cdf, Erlang.CDF(shape, invScale, x), 14); } } }