From cef08e67d17c0f8836e1de332ccafee67ffeb7fc Mon Sep 17 00:00:00 2001 From: Christoph Ruegg Date: Sat, 24 Aug 2013 12:21:16 +0200 Subject: [PATCH] Distributions: adapt Exponential, add InvCDF --- src/Numerics/Distributions/Exponential.cs | 128 ++++++++++++------ .../Continuous/ExponentialTests.cs | 31 ++++- 2 files changed, 115 insertions(+), 44 deletions(-) diff --git a/src/Numerics/Distributions/Exponential.cs b/src/Numerics/Distributions/Exponential.cs index 5280144d..d06740e2 100644 --- a/src/Numerics/Distributions/Exponential.cs +++ b/src/Numerics/Distributions/Exponential.cs @@ -80,16 +80,6 @@ namespace MathNet.Numerics.Distributions return "Exponential(λ = " + _rate + ")"; } - /// - /// Checks whether the parameters of the distribution are valid. - /// - /// The rate (λ) parameter of the distribution. Range: λ ≥ 0. - /// true when the parameters are valid, false otherwise. - static bool IsValidParameterSet(double rate) - { - return rate >= 0.0; - } - /// /// Sets the parameters of the distribution after checking their validity. /// @@ -97,7 +87,7 @@ namespace MathNet.Numerics.Distributions /// When the parameters are out of range. void SetParameters(double rate) { - if (Control.CheckDistributionParameters && !IsValidParameterSet(rate)) + if (rate < 0.0 || Double.IsNaN(rate)) { throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters); } @@ -200,14 +190,10 @@ namespace MathNet.Numerics.Distributions /// /// The location at which to compute the density. /// the density at . + /// public double Density(double x) { - if (x >= 0.0) - { - return _rate*Math.Exp(-_rate*x); - } - - return 0.0; + return x < 0.0 ? 0.0 : _rate*Math.Exp(-_rate*x); } /// @@ -215,6 +201,7 @@ namespace MathNet.Numerics.Distributions /// /// The location at which to compute the log density. /// the log density at . + /// public double DensityLn(double x) { return Math.Log(_rate) - (_rate*x); @@ -225,14 +212,43 @@ 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 (x >= 0.0) + return x < 0.0 ? 0.0 : 1.0 - Math.Exp(-_rate*x); + } + + /// + /// Computes the inverse of the cumulative distribution function (InvCDF) for the distribution + /// at the given probability. This is also known as the quantile or percent point function. + /// + /// The location at which to compute the inverse cumulative density. + /// the inverse cumulative density at . + /// + public double InverseCumulativeDistribution(double p) + { + return p >= 1.0 ? double.PositiveInfinity : -Math.Log(1 - p)/_rate; + } + + /// + /// Draws a random sample from the distribution. + /// + /// A random number from this distribution. + public double Sample() + { + return SampleUnchecked(_random, _rate); + } + + /// + /// Generates a sequence of samples from the Exponential distribution. + /// + /// a sequence of samples from the distribution. + public IEnumerable Samples() + { + while (true) { - return 1.0 - Math.Exp(-_rate*x); + yield return SampleUnchecked(_random, _rate); } - - return 0.0; } /// @@ -249,28 +265,64 @@ namespace MathNet.Numerics.Distributions r = rnd.NextDouble(); } - return -Math.Log(r)/rate; + return -Math.Log(r) / rate; } /// - /// Draws a random sample from the distribution. + /// Computes the probability density of the distribution (PDF) at x, i.e. ∂P(X ≤ x)/∂x. /// - /// A random number from this distribution. - public double Sample() + /// The rate (λ) parameter of the distribution. Range: λ ≥ 0. + /// The location at which to compute the density. + /// the density at . + /// + public static double PDF(double rate, double x) { - return SampleUnchecked(_random, _rate); + if (rate < 0.0) throw new ArgumentOutOfRangeException("rate", Resources.InvalidDistributionParameters); + + return x < 0.0 ? 0.0 : rate*Math.Exp(-rate*x); } /// - /// Generates a sequence of samples from the Exponential 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 rate (λ) parameter of the distribution. Range: λ ≥ 0. + /// The location at which to compute the density. + /// the log density at . + /// + public static double PDFLn(double rate, double x) { - while (true) - { - yield return SampleUnchecked(_random, _rate); - } + if (rate < 0.0) throw new ArgumentOutOfRangeException("rate", Resources.InvalidDistributionParameters); + + return Math.Log(rate) - (rate*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 rate (λ) parameter of the distribution. Range: λ ≥ 0. + /// the cumulative distribution at location . + /// + public static double CDF(double rate, double x) + { + if (rate < 0.0) throw new ArgumentOutOfRangeException("rate", Resources.InvalidDistributionParameters); + + return x < 0.0 ? 0.0 : 1.0 - Math.Exp(-rate*x); + } + + /// + /// Computes the inverse of the cumulative distribution function (InvCDF) for the distribution + /// at the given probability. This is also known as the quantile or percent point function. + /// + /// The location at which to compute the inverse cumulative density. + /// The rate (λ) parameter of the distribution. Range: λ ≥ 0. + /// the inverse cumulative density at . + /// + public static double InvCDF(double rate, double p) + { + if (rate < 0.0) throw new ArgumentOutOfRangeException("rate", Resources.InvalidDistributionParameters); + + return p >= 1.0 ? double.PositiveInfinity : -Math.Log(1 - p)/rate; } /// @@ -281,10 +333,7 @@ namespace MathNet.Numerics.Distributions /// A random number from this distribution. public static double Sample(System.Random rnd, double rate) { - if (Control.CheckDistributionParameters && !IsValidParameterSet(rate)) - { - throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters); - } + if (rate < 0.0) throw new ArgumentOutOfRangeException("rate", Resources.InvalidDistributionParameters); return SampleUnchecked(rnd, rate); } @@ -297,10 +346,7 @@ namespace MathNet.Numerics.Distributions /// a sequence of samples from the distribution. public static IEnumerable Samples(System.Random rnd, double rate) { - if (Control.CheckDistributionParameters && !IsValidParameterSet(rate)) - { - throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters); - } + if (rate < 0.0) throw new ArgumentOutOfRangeException("rate", Resources.InvalidDistributionParameters); while (true) { diff --git a/src/UnitTests/DistributionTests/Continuous/ExponentialTests.cs b/src/UnitTests/DistributionTests/Continuous/ExponentialTests.cs index 5205dcd5..ab55a8a1 100644 --- a/src/UnitTests/DistributionTests/Continuous/ExponentialTests.cs +++ b/src/UnitTests/DistributionTests/Continuous/ExponentialTests.cs @@ -266,11 +266,13 @@ namespace MathNet.Numerics.UnitTests.DistributionTests.Continuous var n = new Exponential(lambda); if (x >= 0) { - Assert.AreEqual(lambda * Math.Exp(-lambda * x), n.Density(x)); + Assert.AreEqual(lambda*Math.Exp(-lambda*x), n.Density(x)); + Assert.AreEqual(lambda*Math.Exp(-lambda*x), Exponential.PDF(lambda, x)); } else { Assert.AreEqual(0.0, n.Density(lambda)); + Assert.AreEqual(0.0, Exponential.PDF(lambda, lambda)); } } @@ -302,7 +304,8 @@ namespace MathNet.Numerics.UnitTests.DistributionTests.Continuous public void ValidateDensityLn(double lambda, double x) { var n = new Exponential(lambda); - Assert.AreEqual(Math.Log(lambda) - (lambda * x), n.DensityLn(x)); + Assert.AreEqual(Math.Log(lambda) - (lambda*x), n.DensityLn(x)); + Assert.AreEqual(Math.Log(lambda) - (lambda*x), Exponential.PDFLn(lambda, x)); } /// @@ -356,12 +359,34 @@ namespace MathNet.Numerics.UnitTests.DistributionTests.Continuous var n = new Exponential(lambda); if (x >= 0.0) { - Assert.AreEqual(1.0 - Math.Exp(-lambda * x), n.CumulativeDistribution(x)); + Assert.AreEqual(1.0 - Math.Exp(-lambda*x), n.CumulativeDistribution(x)); + Assert.AreEqual(1.0 - Math.Exp(-lambda*x), Exponential.CDF(lambda, x)); } else { Assert.AreEqual(0.0, n.CumulativeDistribution(x)); + Assert.AreEqual(0.0, Exponential.CDF(lambda, x)); } } + + /// + /// Validate inverse cumulative distribution. + /// + /// Lambda value. + /// Input X value. + [TestCase(0.1, 0.0)] + [TestCase(1.0, 0.0)] + [TestCase(10.0, 0.0)] + [TestCase(10.0, 0.1)] + [TestCase(1.0, 1.0)] + [TestCase(0.1, Double.PositiveInfinity)] + [TestCase(1.0, Double.PositiveInfinity)] + [TestCase(10.0, Double.PositiveInfinity)] + public void ValidateInverseCumulativeDistribution(double lambda, double x) + { + var n = new Exponential(lambda); + Assert.AreEqual(x, n.InverseCumulativeDistribution(1.0 - Math.Exp(-lambda*x))); + Assert.AreEqual(x, Exponential.InvCDF(lambda, 1.0 - Math.Exp(-lambda*x))); + } } }