From 221da13f427d13819b8eba0b0ff6495ec59b8cfa Mon Sep 17 00:00:00 2001 From: Christoph Ruegg Date: Sat, 24 Aug 2013 09:47:50 +0200 Subject: [PATCH] Distributions: adapt Cauchy, add InvCDF --- src/Numerics/Distributions/Cauchy.cs | 123 +++++++++++++----- .../Continuous/CauchyTests.cs | 39 +++++- 2 files changed, 129 insertions(+), 33 deletions(-) diff --git a/src/Numerics/Distributions/Cauchy.cs b/src/Numerics/Distributions/Cauchy.cs index d0e016dd..bea65ade 100644 --- a/src/Numerics/Distributions/Cauchy.cs +++ b/src/Numerics/Distributions/Cauchy.cs @@ -3,7 +3,9 @@ // http://numerics.mathdotnet.com // http://github.com/mathnet/mathnet-numerics // http://mathnetnumerics.codeplex.com -// Copyright (c) 2009-2010 Math.NET +// +// Copyright (c) 2009-2013 Math.NET +// // Permission is hereby granted, free of charge, to any person // obtaining a copy of this software and associated documentation // files (the "Software"), to deal in the Software without @@ -12,8 +14,10 @@ // copies of the Software, and to permit persons to whom the // Software is furnished to do so, subject to the following // conditions: +// // The above copyright notice and this permission notice shall be // included in all copies or substantial portions of the Software. +// // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, // EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES // OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND @@ -86,17 +90,6 @@ namespace MathNet.Numerics.Distributions return "Cauchy(x0 = " + _location + ", γ = " + _scale + ")"; } - /// - /// Checks whether the parameters of the distribution are valid. - /// - /// The location (x0) of the distribution. - /// The scale (γ) of the distribution. Range: γ > 0. - /// True when the parameters are valid, false otherwise. - static bool IsValidParameterSet(double location, double scale) - { - return scale > 0.0 && !Double.IsNaN(location); - } - /// /// Sets the parameters of the distribution after checking their validity. /// @@ -105,7 +98,7 @@ namespace MathNet.Numerics.Distributions /// When the parameters are out of range. void SetParameters(double location, double scale) { - if (Control.CheckDistributionParameters && !IsValidParameterSet(location, scale)) + if (scale <= 0.0 || Double.IsNaN(location) || Double.IsNaN(scale)) { throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters); } @@ -218,6 +211,7 @@ namespace MathNet.Numerics.Distributions /// /// The location at which to compute the density. /// the density at . + /// public double Density(double x) { return 1.0/(Constants.Pi*_scale*(1.0 + (((x - _location)/_scale)*((x - _location)/_scale)))); @@ -228,6 +222,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(Constants.Pi*_scale*(1.0 + (((x - _location)/_scale)*((x - _location)/_scale)))); @@ -238,22 +233,23 @@ 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/Constants.Pi)*Math.Atan((x - _location)/_scale)) + 0.5; } /// - /// Samples the distribution. + /// 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 location (x0) of the distribution. - /// The scale (γ) of the distribution. Range: γ > 0. - /// a random number from the distribution. - static double SampleUnchecked(System.Random rnd, double location, double scale) + /// The location at which to compute the inverse cumulative density. + /// the inverse cumulative density at . + /// + public double InverseCumulativeDistribution(double p) { - var u = rnd.NextDouble(); - return location + (scale*Math.Tan(Constants.Pi*(u - 0.5))); + return p <= 0.0 ? double.NegativeInfinity : p >= 1.0 ? double.PositiveInfinity + : _location + _scale*Math.Tan((p - 0.5)*Constants.Pi); } /// @@ -277,6 +273,81 @@ namespace MathNet.Numerics.Distributions } } + /// + /// Samples the distribution. + /// + /// The random number generator to use. + /// The location (x0) of the distribution. + /// The scale (γ) of the distribution. Range: γ > 0. + /// a random number from the distribution. + static double SampleUnchecked(System.Random rnd, double location, double scale) + { + var u = rnd.NextDouble(); + return location + (scale*Math.Tan(Constants.Pi*(u - 0.5))); + } + + /// + /// Computes the probability density of the distribution (PDF) at x, i.e. ∂P(X ≤ x)/∂x. + /// + /// The location (x0) of the distribution. + /// The scale (γ) of the distribution. Range: γ > 0. + /// The location at which to compute the density. + /// the density at . + /// + public static double PDF(double location, double scale, double x) + { + if (scale <= 0.0) throw new ArgumentOutOfRangeException("scale", Resources.InvalidDistributionParameters); + + return 1.0/(Constants.Pi*scale*(1.0 + (((x - location)/scale)*((x - location)/scale)))); + } + + /// + /// Computes the log probability density of the distribution (lnPDF) at x, i.e. ln(∂P(X ≤ x)/∂x). + /// + /// The location (x0) of the distribution. + /// The scale (γ) of the distribution. Range: γ > 0. + /// The location at which to compute the density. + /// the log density at . + /// + public static double PDFLn(double location, double scale, double x) + { + if (scale <= 0.0) throw new ArgumentOutOfRangeException("scale", Resources.InvalidDistributionParameters); + + return -Math.Log(Constants.Pi*scale*(1.0 + (((x - location)/scale)*((x - location)/scale)))); + } + + /// + /// 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 location (x0) of the distribution. + /// The scale (γ) of the distribution. Range: γ > 0. + /// the cumulative distribution at location . + /// + public static double CDF(double location, double scale, double x) + { + if (scale <= 0.0) throw new ArgumentOutOfRangeException("scale", Resources.InvalidDistributionParameters); + + return Math.Atan((x - location)/scale)/Constants.Pi + 0.5; + } + + /// + /// 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 location (x0) of the distribution. + /// The scale (γ) of the distribution. Range: γ > 0. + /// the inverse cumulative density at . + /// + public static double InvCDF(double location, double scale, double p) + { + if (scale <= 0.0) throw new ArgumentOutOfRangeException("scale", Resources.InvalidDistributionParameters); + + return p <= 0.0 ? double.NegativeInfinity : p >= 1.0 ? double.PositiveInfinity + : location + scale*Math.Tan((p - 0.5)*Constants.Pi); + } + /// /// Generates a sample from the distribution. /// @@ -286,10 +357,7 @@ namespace MathNet.Numerics.Distributions /// a sample from the distribution. public static double Sample(System.Random rnd, double location, double scale) { - if (Control.CheckDistributionParameters && !IsValidParameterSet(location, scale)) - { - throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters); - } + if (scale <= 0.0) throw new ArgumentOutOfRangeException("scale", Resources.InvalidDistributionParameters); return SampleUnchecked(rnd, location, scale); } @@ -303,10 +371,7 @@ namespace MathNet.Numerics.Distributions /// a sequence of samples from the distribution. public static IEnumerable Samples(System.Random rnd, double location, double scale) { - if (Control.CheckDistributionParameters && !IsValidParameterSet(location, scale)) - { - throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters); - } + if (scale <= 0.0) throw new ArgumentOutOfRangeException("scale", Resources.InvalidDistributionParameters); while (true) { diff --git a/src/UnitTests/DistributionTests/Continuous/CauchyTests.cs b/src/UnitTests/DistributionTests/Continuous/CauchyTests.cs index 0c9c8396..e0f4a0a3 100644 --- a/src/UnitTests/DistributionTests/Continuous/CauchyTests.cs +++ b/src/UnitTests/DistributionTests/Continuous/CauchyTests.cs @@ -3,7 +3,9 @@ // http://numerics.mathdotnet.com // http://github.com/mathnet/mathnet-numerics // http://mathnetnumerics.codeplex.com -// Copyright (c) 2009-2010 Math.NET +// +// Copyright (c) 2009-2013 Math.NET +// // Permission is hereby granted, free of charge, to any person // obtaining a copy of this software and associated documentation // files (the "Software"), to deal in the Software without @@ -12,8 +14,10 @@ // copies of the Software, and to permit persons to whom the // Software is furnished to do so, subject to the following // conditions: +// // The above copyright notice and this permission notice shall be // included in all copies or substantial portions of the Software. +// // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, // EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES // OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND @@ -262,7 +266,9 @@ namespace MathNet.Numerics.UnitTests.DistributionTests.Continuous for (var i = 0; i < 11; i++) { var x = i - 5.0; - Assert.AreEqual(1.0 / ((Constants.Pi * scale) * (1.0 + (((x - location) / scale) * ((x - location) / scale)))), n.Density(x)); + double expected = 1.0 / ((Constants.Pi * scale) * (1.0 + (((x - location) / scale) * ((x - location) / scale)))); + Assert.AreEqual(expected, n.Density(x)); + Assert.AreEqual(expected, Cauchy.PDF(location, scale, x)); } } @@ -283,7 +289,9 @@ namespace MathNet.Numerics.UnitTests.DistributionTests.Continuous for (var i = 0; i < 11; i++) { var x = i - 5.0; - Assert.AreEqual(-Math.Log((Constants.Pi * scale) * (1.0 + (((x - location) / scale) * ((x - location) / scale)))), n.DensityLn(x)); + double expected = -Math.Log((Constants.Pi * scale) * (1.0 + (((x - location) / scale) * ((x - location) / scale)))); + Assert.AreEqual(expected, n.DensityLn(x)); + Assert.AreEqual(expected, Cauchy.PDFLn(location, scale, x)); } } @@ -324,7 +332,30 @@ namespace MathNet.Numerics.UnitTests.DistributionTests.Continuous for (var i = 0; i < 11; i++) { var x = i - 5.0; - Assert.AreEqual(((1.0 / Constants.Pi) * Math.Atan((x - location) / scale)) + 0.5, n.CumulativeDistribution(x)); + double expected = (Math.Atan((x - location)/scale))/Math.PI + 0.5; + Assert.AreEqual(expected, n.CumulativeDistribution(x), 1e-12); + Assert.AreEqual(expected, Cauchy.CDF(location, scale, x), 1e-12); + } + } + + /// + /// Validate inverse cumulative distribution. + /// + /// Location value. + /// Scale value. + [TestCase(0.0, 0.1)] + [TestCase(0.0, 1.0)] + [TestCase(0.0, 10.0)] + [TestCase(-5.0, 100.0)] + public void ValidateInverseCumulativeDistribution(double location, double scale) + { + var n = new Cauchy(location, scale); + for (var i = 0; i < 11; i++) + { + var x = i - 5.0; + double expected = (Math.Atan((x - location)/scale))/Math.PI + 0.5; + Assert.AreEqual(x, n.InverseCumulativeDistribution(expected), 1e-12); + Assert.AreEqual(x, Cauchy.InvCDF(location, scale, expected), 1e-12); } } }