From b8c68e977564e2f1fcde02ddf15dcd26035ad0d4 Mon Sep 17 00:00:00 2001 From: Christoph Ruegg Date: Sat, 24 Aug 2013 16:18:36 +0200 Subject: [PATCH] Distributions: adapt Pareto, add InvCDF --- src/Numerics/Distributions/Pareto.cs | 113 +++++++++++++----- .../Continuous/ParetoTests.cs | 65 +++++----- 2 files changed, 111 insertions(+), 67 deletions(-) diff --git a/src/Numerics/Distributions/Pareto.cs b/src/Numerics/Distributions/Pareto.cs index 2caf58d5..3b546ac1 100644 --- a/src/Numerics/Distributions/Pareto.cs +++ b/src/Numerics/Distributions/Pareto.cs @@ -87,17 +87,6 @@ namespace MathNet.Numerics.Distributions return "Pareto(xm = " + _scale + ", α = " + _shape + ")"; } - /// - /// Checks whether the parameters of the distribution are valid. - /// - /// The scale (xm) of the distribution. Range: xm > 0. - /// The shape (α) of the distribution. Range: α > 0. - /// true when the parameters are valid, false otherwise. - static bool IsValidParameterSet(double scale, double shape) - { - return scale > 0.0 && shape > 0.0; - } - /// /// Sets the parameters of the distribution after checking their validity. /// @@ -106,7 +95,7 @@ namespace MathNet.Numerics.Distributions /// When the parameters are out of range. void SetParameters(double scale, double shape) { - if (Control.CheckDistributionParameters && !IsValidParameterSet(scale, shape)) + if (scale <= 0.0 || shape <= 0.0 || Double.IsNaN(scale) || Double.IsNaN(shape)) { throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters); } @@ -235,6 +224,7 @@ namespace MathNet.Numerics.Distributions /// /// The location at which to compute the density. /// the density at . + /// public double Density(double x) { return _shape*Math.Pow(_scale, _shape)/Math.Pow(x, _shape + 1.0); @@ -245,9 +235,10 @@ 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(Density(x)); + return Math.Log(_shape) + _shape*Math.Log(_scale) - (_shape + 1.0)*Math.Log(x); } /// @@ -255,21 +246,22 @@ namespace MathNet.Numerics.Distributions /// /// The location at which to compute the cumulative distribution function. /// the cumulative distribution at location . + /// public double CumulativeDistribution(double x) { return 1.0 - Math.Pow(_scale/x, _shape); } /// - /// Generates a sample from the Pareto distribution without doing parameter checking. + /// 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 random number generator to use. - /// The scale (xm) of the distribution. Range: xm > 0. - /// The shape (α) of the distribution. Range: α > 0. - /// a random number from the Pareto distribution. - static double SampleUnchecked(System.Random rnd, double scale, double shape) + /// The location at which to compute the inverse cumulative density. + /// the inverse cumulative density at . + /// + public double InverseCumulativeDistribution(double p) { - return scale*Math.Pow(rnd.NextDouble(), -1.0/shape); + return _scale*Math.Pow(1.0 - p, -1.0/_shape); } /// @@ -278,7 +270,7 @@ namespace MathNet.Numerics.Distributions /// A random number from this distribution. public double Sample() { - return SampleUnchecked(_random, _scale, _shape); + return _scale*Math.Pow(_random.NextDouble(), -1.0/_shape); } /// @@ -287,12 +279,74 @@ namespace MathNet.Numerics.Distributions /// a sequence of samples from the distribution. public IEnumerable Samples() { + var power = -1.0/_shape; while (true) { - yield return SampleUnchecked(_random, _scale, _shape); + yield return _scale*Math.Pow(_random.NextDouble(), power); } } + /// + /// Computes the probability density of the distribution (PDF) at x, i.e. ∂P(X ≤ x)/∂x. + /// + /// The scale (xm) of the distribution. Range: xm > 0. + /// The shape (α) of the distribution. Range: α > 0. + /// The location at which to compute the density. + /// the density at . + /// + public static double PDF(double scale, double shape, double x) + { + if (scale <= 0.0 || shape <= 0.0) throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters); + + return shape*Math.Pow(scale, shape)/Math.Pow(x, shape + 1.0); + } + + /// + /// Computes the log probability density of the distribution (lnPDF) at x, i.e. ln(∂P(X ≤ x)/∂x). + /// + /// The scale (xm) of the distribution. Range: xm > 0. + /// The shape (α) of the distribution. Range: α > 0. + /// The location at which to compute the density. + /// the log density at . + /// + public static double PDFLn(double scale, double shape, double x) + { + if (scale <= 0.0 || shape <= 0.0) throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters); + + return Math.Log(shape) + shape*Math.Log(scale) - (shape + 1.0)*Math.Log(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 scale (xm) of the distribution. Range: xm > 0. + /// The shape (α) of the distribution. Range: α > 0. + /// the cumulative distribution at location . + /// + public static double CDF(double scale, double shape, double x) + { + if (scale <= 0.0 || shape <= 0.0) throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters); + + return 1.0 - Math.Pow(scale/x, shape); + } + + /// + /// 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 scale (xm) of the distribution. Range: xm > 0. + /// The shape (α) of the distribution. Range: α > 0. + /// the inverse cumulative density at . + /// + public static double InvCDF(double scale, double shape, double p) + { + if (scale <= 0.0 || shape <= 0.0) throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters); + + return scale*Math.Pow(1.0 - p, -1.0/shape); + } + /// /// Generates a sample from the distribution. /// @@ -302,12 +356,9 @@ namespace MathNet.Numerics.Distributions /// a sample from the distribution. public static double Sample(System.Random rnd, double scale, double shape) { - if (Control.CheckDistributionParameters && !IsValidParameterSet(scale, shape)) - { - throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters); - } + if (scale <= 0.0 || shape <= 0.0) throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters); - return SampleUnchecked(rnd, scale, shape); + return scale*Math.Pow(rnd.NextDouble(), -1.0/shape); } /// @@ -319,14 +370,12 @@ namespace MathNet.Numerics.Distributions /// a sequence of samples from the distribution. public static IEnumerable Samples(System.Random rnd, double scale, double shape) { - if (Control.CheckDistributionParameters && !IsValidParameterSet(scale, shape)) - { - throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters); - } + if (scale <= 0.0 || shape <= 0.0) throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters); + var power = -1.0 / shape; while (true) { - yield return SampleUnchecked(rnd, scale, shape); + yield return scale*Math.Pow(rnd.NextDouble(), power); } } } diff --git a/src/UnitTests/DistributionTests/Continuous/ParetoTests.cs b/src/UnitTests/DistributionTests/Continuous/ParetoTests.cs index 8f554a45..d4f17cff 100644 --- a/src/UnitTests/DistributionTests/Continuous/ParetoTests.cs +++ b/src/UnitTests/DistributionTests/Continuous/ParetoTests.cs @@ -282,40 +282,53 @@ namespace MathNet.Numerics.UnitTests.DistributionTests.Continuous Assert.AreEqual(Double.PositiveInfinity, n.Maximum); } - /// - /// Validate density. - /// - /// Scale value. - /// Shape value. - /// Input X value. + [TestCase(1, 1, 1, 1)] + [TestCase(1, 1, 1.5, 4/9.0)] + [TestCase(1, 1, 5, 1/25.0)] + [TestCase(1, 1, 50, 1/2500.0)] + [TestCase(1, 4, 1, 4)] + [TestCase(1, 4, 1.5, 128/243.0)] + [TestCase(1, 4, 50, 1/78125000.0)] + [TestCase(3, 2, 3, 2/3.0)] + [TestCase(3, 2, 5, 18/125.0)] + [TestCase(25, 100, 50, 1.5777218104420236e-30)] + [TestCase(100, 25, 150, 6.6003546737276816e-6)] + public void ValidateDensity(double scale, double shape, double x, double expected) + { + var dist = new Pareto(scale, shape); + + Assert.AreEqual(expected, dist.Density(x), 1e-12); + Assert.AreEqual(expected, Pareto.PDF(scale, shape, x), 1e-12); + + Assert.AreEqual(Math.Log(expected), dist.DensityLn(x), 1e-12); + Assert.AreEqual(Math.Log(expected), Pareto.PDFLn(scale, shape, x), 1e-12); + } + [TestCase(0.1, 0.1, 0.1)] [TestCase(1.0, 1.0, 1.0)] [TestCase(5.0, 5.0, 2.0)] [TestCase(7.0, 7.0, 10.0)] [TestCase(10.0, 10.0, 12.0)] [TestCase(Double.PositiveInfinity, Double.PositiveInfinity, Double.PositiveInfinity)] - public void ValidateDensity(double scale, double shape, double x) + public void ValidateCumulativeDistribution(double scale, double shape, double x) { var n = new Pareto(scale, shape); - Assert.AreEqual(shape * Math.Pow(scale, shape) / Math.Pow(x, shape + 1.0), n.Density(x)); + double expected = 1.0 - Math.Pow(scale/x, shape); + Assert.AreEqual(expected, n.CumulativeDistribution(x)); + Assert.AreEqual(expected, Pareto.CDF(scale, shape, x)); } - /// - /// Validate density log. - /// - /// Scale value. - /// Shape value. - /// Input X value. [TestCase(0.1, 0.1, 0.1)] [TestCase(1.0, 1.0, 1.0)] [TestCase(5.0, 5.0, 2.0)] [TestCase(7.0, 7.0, 10.0)] [TestCase(10.0, 10.0, 12.0)] - [TestCase(Double.PositiveInfinity, Double.PositiveInfinity, Double.PositiveInfinity)] - public void ValidateDensityLn(double scale, double shape, double x) + public void ValidateInverseCumulativeDistribution(double scale, double shape, double x) { var n = new Pareto(scale, shape); - Assert.AreEqual(Math.Log(n.Density(x)), n.DensityLn(x)); + double cdf = 1.0 - Math.Pow(scale / x, shape); + Assert.AreEqual(x, n.InverseCumulativeDistribution(cdf), 1e-12); + Assert.AreEqual(x, Pareto.InvCDF(scale, shape, cdf), 1e-12); } /// @@ -338,23 +351,5 @@ namespace MathNet.Numerics.UnitTests.DistributionTests.Continuous var ied = n.Samples(); ied.Take(5).ToArray(); } - - /// - /// Validate cumulative distribution. - /// - /// Scale value. - /// Shape value. - /// Input X value. - [TestCase(0.1, 0.1, 0.1)] - [TestCase(1.0, 1.0, 1.0)] - [TestCase(5.0, 5.0, 2.0)] - [TestCase(7.0, 7.0, 10.0)] - [TestCase(10.0, 10.0, 12.0)] - [TestCase(Double.PositiveInfinity, Double.PositiveInfinity, Double.PositiveInfinity)] - public void ValidateCumulativeDistribution(double scale, double shape, double x) - { - var n = new Pareto(scale, shape); - Assert.AreEqual(1.0 - Math.Pow(scale / x, shape), n.CumulativeDistribution(x)); - } } }