Distributions: adapt Pareto, add InvCDF

src/Numerics/Distributions/Pareto.cs

@@ -87,17 +87,6 @@ namespace MathNet.Numerics.Distributions
             return "Pareto(xm = " + _scale + ", α = " + _shape + ")";
         }
 
-        /// <summary>
-        /// Checks whether the parameters of the distribution are valid. 
-        /// </summary>
-        /// <param name="scale">The scale (xm) of the distribution. Range: xm > 0.</param>
-        /// <param name="shape">The shape (α) of the distribution. Range: α > 0.</param>
-        /// <returns><c>true</c> when the parameters are valid, <c>false</c> otherwise.</returns>
-        static bool IsValidParameterSet(double scale, double shape)
-        {
-            return scale > 0.0 && shape > 0.0;
-        }
-
         /// <summary>
         /// Sets the parameters of the distribution after checking their validity.
         /// </summary>
@@ -106,7 +95,7 @@ namespace MathNet.Numerics.Distributions
         /// <exception cref="ArgumentOutOfRangeException">When the parameters are out of range.</exception>
         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
         /// </summary>
         /// <param name="x">The location at which to compute the density.</param>
         /// <returns>the density at <paramref name="x"/>.</returns>
+        /// <seealso cref="PDF"/>
         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
         /// </summary>
         /// <param name="x">The location at which to compute the log density.</param>
         /// <returns>the log density at <paramref name="x"/>.</returns>
+        /// <seealso cref="PDFLn"/>
         public double DensityLn(double x)
         {
-            return Math.Log(Density(x));
+            return Math.Log(_shape) + _shape*Math.Log(_scale) - (_shape + 1.0)*Math.Log(x);
         }
 
         /// <summary>
@@ -255,21 +246,22 @@ namespace MathNet.Numerics.Distributions
         /// </summary>
         /// <param name="x">The location at which to compute the cumulative distribution function.</param>
         /// <returns>the cumulative distribution at location <paramref name="x"/>.</returns>
+        /// <seealso cref="CDF"/>
         public double CumulativeDistribution(double x)
         {
             return 1.0 - Math.Pow(_scale/x, _shape);
         }
 
         /// <summary>
-        /// 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.
         /// </summary>
-        /// <param name="rnd">The random number generator to use.</param>
+        /// <param name="p">The location at which to compute the inverse cumulative density.</param>
-        /// <param name="scale">The scale (xm) of the distribution. Range: xm > 0.</param>
+        /// <returns>the inverse cumulative density at <paramref name="p"/>.</returns>
-        /// <param name="shape">The shape (α) of the distribution. Range: α > 0.</param>
+        /// <seealso cref="InvCDF"/>
-        /// <returns>a random number from the Pareto distribution.</returns>
+        public double InverseCumulativeDistribution(double p)
-        static double SampleUnchecked(System.Random rnd, double scale, double shape)
         {
-            return scale*Math.Pow(rnd.NextDouble(), -1.0/shape);
+            return _scale*Math.Pow(1.0 - p, -1.0/_shape);
         }
 
         /// <summary>
@@ -278,7 +270,7 @@ namespace MathNet.Numerics.Distributions
         /// <returns>A random number from this distribution.</returns>
         public double Sample()
         {
-            return SampleUnchecked(_random, _scale, _shape);
+            return _scale*Math.Pow(_random.NextDouble(), -1.0/_shape);
         }
 
         /// <summary>
@@ -287,12 +279,74 @@ namespace MathNet.Numerics.Distributions
         /// <returns>a sequence of samples from the distribution.</returns>
         public IEnumerable<double> Samples()
         {
+            var power = -1.0/_shape;
             while (true)
             {
-                yield return SampleUnchecked(_random, _scale, _shape);
+                yield return _scale*Math.Pow(_random.NextDouble(), power);
             }
         }
 
+        /// <summary>
+        /// Computes the probability density of the distribution (PDF) at x, i.e. ∂P(X ≤ x)/∂x.
+        /// </summary>
+        /// <param name="scale">The scale (xm) of the distribution. Range: xm > 0.</param>
+        /// <param name="shape">The shape (α) of the distribution. Range: α > 0.</param>
+        /// <param name="x">The location at which to compute the density.</param>
+        /// <returns>the density at <paramref name="x"/>.</returns>
+        /// <seealso cref="Density"/>
+        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);
+        }
+
+        /// <summary>
+        /// Computes the log probability density of the distribution (lnPDF) at x, i.e. ln(∂P(X ≤ x)/∂x).
+        /// </summary>
+        /// <param name="scale">The scale (xm) of the distribution. Range: xm > 0.</param>
+        /// <param name="shape">The shape (α) of the distribution. Range: α > 0.</param>
+        /// <param name="x">The location at which to compute the density.</param>
+        /// <returns>the log density at <paramref name="x"/>.</returns>
+        /// <seealso cref="DensityLn"/>
+        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);
+        }
+
+        /// <summary>
+        /// Computes the cumulative distribution (CDF) of the distribution at x, i.e. P(X ≤ x).
+        /// </summary>
+        /// <param name="x">The location at which to compute the cumulative distribution function.</param>
+        /// <param name="scale">The scale (xm) of the distribution. Range: xm > 0.</param>
+        /// <param name="shape">The shape (α) of the distribution. Range: α > 0.</param>
+        /// <returns>the cumulative distribution at location <paramref name="x"/>.</returns>
+        /// <seealso cref="CumulativeDistribution"/>
+        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);
+        }
+
+        /// <summary>
+        /// 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.
+        /// </summary>
+        /// <param name="p">The location at which to compute the inverse cumulative density.</param>
+        /// <param name="scale">The scale (xm) of the distribution. Range: xm > 0.</param>
+        /// <param name="shape">The shape (α) of the distribution. Range: α > 0.</param>
+        /// <returns>the inverse cumulative density at <paramref name="p"/>.</returns>
+        /// <seealso cref="InverseCumulativeDistribution"/>
+        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);
+        }
+
         /// <summary>
         /// Generates a sample from the distribution.
         /// </summary>
@@ -302,12 +356,9 @@ namespace MathNet.Numerics.Distributions
         /// <returns>a sample from the distribution.</returns>
         public static double Sample(System.Random rnd, double scale, double shape)
         {
-            if (Control.CheckDistributionParameters && !IsValidParameterSet(scale, shape))
+            if (scale <= 0.0 || shape <= 0.0) throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters);
-            {
-                throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters);
-            }
 
-            return SampleUnchecked(rnd, scale, shape);
+            return scale*Math.Pow(rnd.NextDouble(), -1.0/shape);
         }
 
         /// <summary>
@@ -319,14 +370,12 @@ namespace MathNet.Numerics.Distributions
         /// <returns>a sequence of samples from the distribution.</returns>
         public static IEnumerable<double> Samples(System.Random rnd, double scale, double shape)
         {
-            if (Control.CheckDistributionParameters && !IsValidParameterSet(scale, shape))
+            if (scale <= 0.0 || shape <= 0.0) throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters);
-            {
-                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);
             }
         }
     }

src/UnitTests/DistributionTests/Continuous/ParetoTests.cs

@@ -282,40 +282,53 @@ namespace MathNet.Numerics.UnitTests.DistributionTests.Continuous
             Assert.AreEqual(Double.PositiveInfinity, n.Maximum);
         }
 
-        /// <summary>
+        [TestCase(1, 1, 1, 1)]
-        /// Validate density.
+        [TestCase(1, 1, 1.5, 4/9.0)]
-        /// </summary>
+        [TestCase(1, 1, 5, 1/25.0)]
-        /// <param name="scale">Scale value.</param>
+        [TestCase(1, 1, 50, 1/2500.0)]
-        /// <param name="shape">Shape value.</param>
+        [TestCase(1, 4, 1, 4)]
-        /// <param name="x">Input X value.</param>
+        [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));
         }
 
-        /// <summary>
-        /// Validate density log.
-        /// </summary>
-        /// <param name="scale">Scale value.</param>
-        /// <param name="shape">Shape value.</param>
-        /// <param name="x">Input X value.</param>
         [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 ValidateInverseCumulativeDistribution(double scale, double shape, double x)
-        public void ValidateDensityLn(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);
         }
 
         /// <summary>
@@ -338,23 +351,5 @@ namespace MathNet.Numerics.UnitTests.DistributionTests.Continuous
             var ied = n.Samples();
             ied.Take(5).ToArray();
         }
-
-        /// <summary>
-        /// Validate cumulative distribution.
-        /// </summary>
-        /// <param name="scale">Scale value.</param>
-        /// <param name="shape">Shape value.</param>
-        /// <param name="x">Input X value.</param>
-        [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));
-        }
     }
 }