From 9028337af9cfcbc41f4212df700cacdf5a007086 Mon Sep 17 00:00:00 2001 From: Philip Gruber <55662520+PGsStudio@users.noreply.github.com> Date: Sun, 6 Oct 2019 21:38:33 +0200 Subject: [PATCH] Added new distributions and corresponding tests - New distributions are - Three parameter Burr (Type XII) - Inverse Gaussian - Truncated Pareto - Unit tests added for all public methods - New distributions added to CommonDistributionTests --- .../CommonDistributionTests.cs | 3 + .../DistributionTests/Continuous/BurrTests.cs | 426 ++++++++++++++++ .../Continuous/InverseGaussianTests.cs | 419 ++++++++++++++++ .../Continuous/TruncatedParetoTests.cs | 446 +++++++++++++++++ src/Numerics/Distributions/Burr.cs | 429 ++++++++++++++++ src/Numerics/Distributions/InverseGaussian.cs | 449 +++++++++++++++++ src/Numerics/Distributions/TruncatedPareto.cs | 466 ++++++++++++++++++ src/Numerics/Precision.Comparison.cs | 9 + src/Numerics/Properties/Resources.Designer.cs | 11 + src/Numerics/Properties/Resources.resx | 3 + 10 files changed, 2661 insertions(+) create mode 100644 src/Numerics.Tests/DistributionTests/Continuous/BurrTests.cs create mode 100644 src/Numerics.Tests/DistributionTests/Continuous/InverseGaussianTests.cs create mode 100644 src/Numerics.Tests/DistributionTests/Continuous/TruncatedParetoTests.cs create mode 100644 src/Numerics/Distributions/Burr.cs create mode 100644 src/Numerics/Distributions/InverseGaussian.cs create mode 100644 src/Numerics/Distributions/TruncatedPareto.cs diff --git a/src/Numerics.Tests/DistributionTests/CommonDistributionTests.cs b/src/Numerics.Tests/DistributionTests/CommonDistributionTests.cs index 7b0c326b..c2fce2c4 100644 --- a/src/Numerics.Tests/DistributionTests/CommonDistributionTests.cs +++ b/src/Numerics.Tests/DistributionTests/CommonDistributionTests.cs @@ -67,6 +67,7 @@ namespace MathNet.Numerics.UnitTests.DistributionTests { new Beta(1.0, 1.0), new BetaScaled(1.0, 1.5, 0.5, 2.0), + new Burr(2.0, 3.0, 1.0), new Cauchy(1.0, 1.0), new Chi(3.0), new ChiSquared(3.0), @@ -76,6 +77,7 @@ namespace MathNet.Numerics.UnitTests.DistributionTests new FisherSnedecor(0.3, 0.4), new Gamma(1.0, 1.0), new InverseGamma(1.0, 1.0), + new InverseGaussian(1.0, 3.0), new Laplace(1.0, 0.5), new LogNormal(1.0, 1.0), new Normal(0.0, 1.0), @@ -84,6 +86,7 @@ namespace MathNet.Numerics.UnitTests.DistributionTests new Stable(0.5, 1.0, 0.5, 1.0), new StudentT(0.0, 1.0, 5.0), new Triangular(0, 1, 0.7), + new TruncatedPareto(1.0, 3.0, 1000), new Weibull(1.0, 1.0), }; diff --git a/src/Numerics.Tests/DistributionTests/Continuous/BurrTests.cs b/src/Numerics.Tests/DistributionTests/Continuous/BurrTests.cs new file mode 100644 index 00000000..d5e7baa3 --- /dev/null +++ b/src/Numerics.Tests/DistributionTests/Continuous/BurrTests.cs @@ -0,0 +1,426 @@ +// +// Math.NET Numerics, part of the Math.NET Project +// http://numerics.mathdotnet.com +// http://github.com/mathnet/mathnet-numerics +// +// Copyright (c) 2009-2019 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 +// restriction, including without limitation the rights to use, +// copy, modify, merge, publish, distribute, sublicense, and/or sell +// 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 +// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT +// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, +// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING +// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR +// OTHER DEALINGS IN THE SOFTWARE. +// + +using MathNet.Numerics.Distributions; +using NUnit.Framework; +using System; +using System.Linq; + +namespace MathNet.Numerics.UnitTests.DistributionTests.Continuous +{ + /// + /// Burr distribution tests. + /// + [TestFixture, Category("Distributions")] + public class BurrTests + { + private const int highPrecision = 12; + private const int lowPrecision = 8; + + /// + /// Can create Burr. + /// + /// a parameter. + /// c parameter. + /// k parameter. + [TestCase(1.0, 1.0, 0.5)] + [TestCase(10.0, 0.1, 0.2)] + [TestCase(5.0, 1.0, 3)] + [TestCase(2.0, 10.0, 2.0)] + [TestCase(10.0, 100.0, 1.0)] + public void CanCreateBurr(double a, double c, double k) + { + var n = new Burr(a, c, k); + Assert.AreEqual(a, n.a); + Assert.AreEqual(c, n.c); + Assert.AreEqual(k, n.k); + } + + /// + /// Can create Burr with random source. + /// + [Test] + public void CanCreateBurrWithRandomSource() + { + var randomSource = new Numerics.Random.MersenneTwister(100); + var n = new Burr(1.0, 1.0, 0.5, randomSource); + Assert.AreEqual(randomSource, n.RandomSource); + } + + /// + /// Burr create fails with bad parameters. + /// + /// a parameter. + /// c parameter. + /// k parameter. + [TestCase(5.0, 1.0, -1.0)] + [TestCase(-1.0, 1.0, 2.05)] + [TestCase(1.0, -1.0, 2.05)] + [TestCase(Double.NaN, Double.NaN, Double.NaN)] + [TestCase(Double.NaN, 1.0, 2.4)] + [TestCase(1.0, 1.0, double.PositiveInfinity)] + [TestCase(1.0, -1.0, Double.NegativeInfinity)] + public void BurrCreateFailsWithBadParameters(double a, double c, double k) + { + Assert.That(() => new Burr(a, c, k), Throws.ArgumentException); + } + + /// + /// Validate IsValidParameterSet. + /// + /// a parameter. + /// c parameter. + /// k parameter. + /// exoected validity of paremeter set. + [TestCase(1.0, 1.0, 0.5, true)] + [TestCase(10.0, 0.1, 0.2, true)] + [TestCase(5.0, 1.0, -1.0, false)] + [TestCase(1.0, 1.0, double.PositiveInfinity, false)] + public void ValidateIsValidParameterSet(double a, double c, double k, bool validity) + { + Assert.AreEqual(Burr.IsValidParameterSet(a, c, k), validity); + } + + /// + /// Validate ToString. + /// + [Test] + public void ValidateToString() + { + var n = new Burr(1d, 2d, 3d); + Assert.AreEqual("Burr(a = 1, c = 2, k = 3)", n.ToString()); + } + + /// + /// Validate mode. + /// + /// a parameter. + /// c parameter. + /// k parameter. + /// Expected value. + [TestCase(2.000000, 2.00000, 1.0, 1.154700538379251)] + [TestCase(1.000000, 1.00000, 0.5, 0.0)] + [TestCase(3.000000, 4.00000, 1.0, 2.640335210380180)] + public void ValidateMode(double a, double c, double k, double mode) + { + var n = new Burr(a, c, k); + AssertHelpers.AlmostEqualRelative(mode, n.Mode, highPrecision); + } + + /// + /// Validate median. + /// + /// a parameter. + /// c parameter. + /// k parameter. + /// Expected value. + [TestCase(1.000000, 1.00000, 0.5, 3.0)] + [TestCase(1.000000, 0.50000, 0.5, 9.0)] + [TestCase(1.000000, 1.00000, 5.0, 0.148698354997035)] + public void ValidateMedian(double a, double c, double k, double median) + { + var n = new Burr(a, c, k); + AssertHelpers.AlmostEqualRelative(median, n.Median, highPrecision); + } + + /// + /// Validate mean. + /// + /// a parameter. + /// c parameter. + /// k parameter. + /// Expected value. + [TestCase(1.000000, 1.00000, 1.5, 2)] + [TestCase(4.000000, 5.00000, 0.5, 6.198785110989412)] + [TestCase(4.000000, 5.00000, 5.0, 2.729694550490384)] + public void ValidateMean(double a, double c, double k, double mean) + { + var n = new Burr(a, c, k); + AssertHelpers.AlmostEqualRelative(mean, n.Mean, highPrecision); + } + + /// + /// Validate variance. + /// + /// a parameter. + /// c parameter. + /// k parameter. + /// Expected value. + [TestCase(4.0, 5.0, 2.0, 0.983559126843161)] + [TestCase(2.0, 3.5, 4.0, 0.207489170174404)] + public void ValidateVariance(double a, double c, double k, double variance) + { + var n = new Burr(a, c, k); + AssertHelpers.AlmostEqualRelative(variance, n.Variance, highPrecision); + } + + /// + /// Validate standard deviation. + /// + /// a parameter. + /// c parameter. + /// k parameter. + /// Expected value. + [TestCase(4.0, 5.0, 2.0, 0.991745494995143)] + [TestCase(2.0, 3.5, 4.0, 0.455509791524182)] + public void ValidateStandardDeviation(double a, double c, double k, double std) + { + var n = new Burr(a, c, k); + AssertHelpers.AlmostEqualRelative(std, n.StdDev, highPrecision); + } + + /// + /// Validate minimum. + /// + [Test] + public void ValidateMinimum() + { + var n = new Burr(1.0, 2.0, 2.0); + Assert.AreEqual(0.0, n.Minimum); + } + + /// + /// Validate maximum. + /// + [Test] + public void ValidateMaximum() + { + var n = new Burr(1.0, 2.0, 1.0); + Assert.AreEqual(Double.PositiveInfinity, n.Maximum); + } + + /// + /// Validate entropy throws not supported exception + /// + /// a parameter. + /// c parameter. + /// k parameter. + [TestCase(1.0, 1.0, 0.5)] + public void ValidateEntropyFailsWithNotSupported(double a, double c, double k) + { + var n = new Burr(a, c, k); + Assert.Throws(() => { var x = n.Entropy; }); + } + + /// + /// Validate GetMoments. + /// + /// a parameter. + /// c parameter. + /// k parameter. + /// order of the moment. + /// Expected value. + [TestCase(1.000000, 1.00000, 1.5, 1, 2)] + [TestCase(4.000000, 5.00000, 0.5, 1, 6.198785110989412)] + [TestCase(4.000000, 5.00000, 5.0, 1, 2.729694550490384)] + [TestCase(4.0, 5.0, 2.0, 3, 50.738165747621750)] + [TestCase(2.0, 3.5, 4.0, 3, 2.895046685294824)] + public void ValidateMoments(double a, double c, double k, int order, double value) + { + var n = new Burr(a, c, k); + AssertHelpers.AlmostEqualRelative(value, n.GetMoment(order), lowPrecision); + } + + /// + /// Validate GetMoments throws when given bad parameters + /// + /// a parameter. + /// c parameter. + /// k parameter. + /// the order of the moment. + [TestCase(1.0, 1.0, 0.5, 1)] + [TestCase(10.0, 1.0, 1.5, 2)] + [TestCase(5.0, 2, 1.4, 3)] + public void ValidateGetMomentsFailsWithBadParameters(double a, double c, double k, double order) + { + var n = new Burr(a, c, k); + Assert.That(() => n.GetMoment(order), Throws.ArgumentException); + } + + /// + /// Validate skewness. + /// + /// a parameter. + /// c parameter. + /// k parameter. + /// Expected value. + [TestCase(4.0, 5.0, 2.0, 0.635277200891842)] + [TestCase(2.0, 3.5, 4.0, 0.483073007212360)] + public void ValidateSkewness(double a, double c, double k, double skewness) + { + var n = new Burr(a, c, k); + AssertHelpers.AlmostEqualRelative(skewness, n.Skewness, lowPrecision); + } + + /// + /// Validate density. + /// + /// a parameter. + /// c parameter. + /// k parameter. + /// Input X value. + /// Expected value. + [TestCase(1.00000, 1.000000, 0.5, 0.100000, 0.433392086020724)] + [TestCase(1.000000, 1.000000, 5.0, 0.500000, 0.438957475994512)] + public void ValidateDensity(double a, double c, double k, double x, double p) + { + var n = new Burr(a, c, k); + AssertHelpers.AlmostEqualRelative(p, n.Density(x), highPrecision); + AssertHelpers.AlmostEqualRelative(p, Burr.PDF(a, c, k, x), highPrecision); + } + + /// + /// Validate density log. + /// + /// a parameter. + /// c parameter. + /// k parameter. + /// Input X value. + /// Expected value. + [TestCase(1.000000, 1.000000, 0.5, 0.500000, -1.301344842722192)] + [TestCase(1.000000, 1.000000, 5.0, 0.500000, -0.823352736214886)] + [TestCase(4.000000, 5.000000, 1.0, 5.000000, -1.682583872895781)] + public void ValidateDensityLn(double a, double c, double k, double x, double p) + { + var n = new Burr(a, c, k); + AssertHelpers.AlmostEqualRelative(p, n.DensityLn(x), highPrecision); + AssertHelpers.AlmostEqualRelative(p, Burr.PDFLn(a, c, k, x), highPrecision); + } + + /// + /// Validate cumulative distribution. + /// + /// a parameter. + /// c parameter. + /// k parameter. + /// The location at which to compute the density. + /// Expected value. + [TestCase(1.000000, 1.000000, 0.5, 0.500000, 0.183503419072274)] + [TestCase(4.000000, 5.000000, 0.5, 0.500000, 0.00001525843981653452)] + [TestCase(4.000000, 5.000000, 0.5, 5.00000, 0.503203804978536)] + [TestCase(4.000000, 5.000000, 5.0, 5.00000, 0.999084238019496)] + public void ValidateCumulativeDistribution(double a, double c, double k, double x, double f) + { + var n = new Burr(a, c, k); + AssertHelpers.AlmostEqualRelative(f, n.CumulativeDistribution(x), lowPrecision); + AssertHelpers.AlmostEqualRelative(f, Burr.CDF(a, c, k, x), lowPrecision); + } + + /// + /// Can sample static. + /// + [Test] + public void CanSampleStatic() + { + Burr.Sample(new Numerics.Random.MersenneTwister(100), 1.0, 2.0, 2.0); + } + + /// + /// Can fill array with samples static. + /// + [Test] + public void CanFillSampleArrayStatic() + { + double[] samples = new double[100]; + Burr.Samples(new Numerics.Random.MersenneTwister(100), samples, 1.0, 2.0, 2.0); + Assert.IsTrue(!samples.Any(x => x == 0)); + } + + /// + /// Can sample sequence static. + /// + [Test] + public void CanSampleSequenceStatic() + { + var ied = Burr.Samples(new Numerics.Random.MersenneTwister(100), 1.0, 2.0, 2.0); + GC.KeepAlive(ied.Take(5).ToArray()); + } + + /// + /// Fail sample static with bad parameters. + /// + [Test] + public void FailSampleStatic() + { + Assert.That(() => { var d = Burr.Sample(new Numerics.Random.MersenneTwister(100), 1.0, -1.0, 2.0); }, Throws.ArgumentException); + } + + /// + /// Fail filling array with samples with bad parameters static. + /// + [Test] + public void FailFillingSampleArrayStatic() + { + double[] samples = new double[100]; + Assert.That(() => { Burr.Samples(new Numerics.Random.MersenneTwister(100), samples, -1.0, 2.0, 2.0); }, Throws.ArgumentException); + } + + /// + /// Fail sample sequence static with bad parameters. + /// + [Test] + public void FailSampleSequenceStatic() + { + Assert.That(() => { var ied = Burr.Samples(new Numerics.Random.MersenneTwister(100), 1.0, -1.0, 2.0).First(); }, Throws.ArgumentException); + } + + /// + /// Can sample. + /// + [Test] + public void CanSample() + { + var n = new Burr(1.0, 2.0, 2.0); + n.Sample(); + } + + /// + /// Can fill array with samples. + /// + [Test] + public void CanFillSampleArray() + { + double[] samples = new double[100]; + var n = new Burr(1.0, 2.0, 2.0, new Numerics.Random.MersenneTwister(100)); + n.Samples(samples); + Assert.IsTrue(!samples.Any(x => x == 0)); + } + + /// + /// Can sample sequence. + /// + [Test] + public void CanSampleSequence() + { + var n = new Burr(1.0, 2.0, 2.0); + var ied = n.Samples(); + GC.KeepAlive(ied.Take(5).ToArray()); + } + } +} diff --git a/src/Numerics.Tests/DistributionTests/Continuous/InverseGaussianTests.cs b/src/Numerics.Tests/DistributionTests/Continuous/InverseGaussianTests.cs new file mode 100644 index 00000000..5c13e79c --- /dev/null +++ b/src/Numerics.Tests/DistributionTests/Continuous/InverseGaussianTests.cs @@ -0,0 +1,419 @@ +// +// Math.NET Numerics, part of the Math.NET Project +// http://numerics.mathdotnet.com +// http://github.com/mathnet/mathnet-numerics +// +// Copyright (c) 2009-2019 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 +// restriction, including without limitation the rights to use, +// copy, modify, merge, publish, distribute, sublicense, and/or sell +// 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 +// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT +// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, +// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING +// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR +// OTHER DEALINGS IN THE SOFTWARE. +// + +using MathNet.Numerics.Distributions; +using NUnit.Framework; +using System; +using System.Linq; + +namespace MathNet.Numerics.UnitTests.DistributionTests.Continuous +{ + /// + /// InverseGaussian distribution tests. + /// + [TestFixture, Category("Distributions")] + public class InverseGaussianTests + { + private const int precision = 12; + + /// + /// Can create Inverse Gaussian. + /// + /// Mu value. + /// Lambda parameter. + [TestCase(0.1, 0.1)] + [TestCase(1.0, 1.0)] + [TestCase(10.0, 10.0)] + public void CanCreateInverseGaussian(double mu, double lambda) + { + var n = new InverseGaussian(mu, lambda); + Assert.AreEqual(mu, n.Mu); + Assert.AreEqual(lambda, n.Lambda); + } + + /// + /// Can create Inverse Gaussian with random source. + /// + [Test] + public void CanCreateInverseGaussianWithRandomSource() + { + var randomSource = new Numerics.Random.MersenneTwister(100); + var n = new InverseGaussian(1.0, 1.0, randomSource); + Assert.AreEqual(randomSource, n.RandomSource); + } + + /// + /// InverseGaussian create fails with bad parameters. + /// + /// Mu parameter. + /// Lambda parameter. + [TestCase(Double.NaN, 1.0)] + [TestCase(1.0, Double.NaN)] + [TestCase(Double.NaN, Double.NaN)] + [TestCase(-1.0, -1.0)] + public void InverseGaussianCreateFailsWithBadParameters(double mu, double lambda) + { + Assert.That(() => new LogNormal(mu, lambda), Throws.ArgumentException); + } + + /// + /// Validate IsValidParameterSet. + /// + /// Mu parameter. + /// Lambda parameter. + /// exoected validity of paremeter set. + [TestCase(1.0, 1.0, true)] + [TestCase(10.0, 10.0, true)] + [TestCase(-1.0, -1.0, false)] + public void ValidateIsValidParameterSet(double mu, double lambda, bool validity) + { + Assert.AreEqual(InverseGaussian.IsValidParameterSet(mu, lambda), validity); + } + + /// + /// Validate ToString. + /// + [Test] + public void ValidateToString() + { + var n = new InverseGaussian(1d, 2d); + Assert.AreEqual("InverseGaussian(μ = 1, λ = 2)", n.ToString()); + } + + /// + /// Validate mode. + /// + /// Mu parameter. + /// Lambda parameter. + /// Expected value. + [TestCase(0.100000, 0.100000, 0.030277563773199)] + [TestCase(1.500000, 5.500000, 1.007026953273370)] + [TestCase(2.500000, 1.500000, 0.481456008918130)] + [TestCase(5.500000, 2.500000, 0.815033614523337)] + [TestCase(5.500000, 5.500000, 1.665266007525970)] + public void ValidateMode(double mu, double lambda, double mode) + { + var n = new InverseGaussian(mu, lambda); + AssertHelpers.AlmostEqualRelative(mode, n.Mode, precision); + } + + /// + /// Validate median. + /// + /// Mu parameter. + /// Lambda parameter. + /// Expected value. + [TestCase(0.100000, 0.100000, 0.067584130569524)] + [TestCase(1.500000, 0.100000, 0.190369896722453)] + [TestCase(2.500000, 0.100000, 0.201110669500521)] + [TestCase(5.500000, 5.500000, 3.717127181323815)] + public void ValidateMedian(double mu, double lambda, double median) + { + var n = new InverseGaussian(mu, lambda); + AssertHelpers.AlmostEqualRelative(median, n.Median, precision); + } + + /// + /// Validate mean. + /// + /// Mu parameter. + /// Lambda parameter. + /// Expected value. + [TestCase(0.100000, 5.500000, 0.1)] + [TestCase(1.500000, 0.100000, 1.5)] + [TestCase(2.500000, 5.500000, 2.5)] + [TestCase(5.500000, 0.100000, 5.5)] + public void ValidateMean(double mu, double lambda, double mean) + { + var n = new InverseGaussian(mu, lambda); + AssertHelpers.AlmostEqualRelative(mean, n.Mean, precision); + } + + /// + /// Validate variance. + /// + /// Mu parameter. + /// Lambda parameter. + /// Expected value. + [TestCase(0.100000, 5.500000, 1.818181818181819e-04)] + [TestCase(1.500000, 0.100000, 33.75)] + public void ValidateVariance(double mu, double lambda, double variance) + { + var n = new InverseGaussian(mu, lambda); + AssertHelpers.AlmostEqualRelative(variance, n.Variance, precision); + } + + /// + /// Validate standard deviation. + /// + /// Mu parameter. + /// Lambda parameter. + /// Expected value. + [TestCase(2.500000, 5.500000, 1.685499656158105)] + [TestCase(5.500000, 0.100000, 40.789091679026143)] + public void ValidateStandardDeviation(double mu, double lambda, double std) + { + var n = new InverseGaussian(mu, lambda); + AssertHelpers.AlmostEqualRelative(std, n.StdDev, precision); + } + + /// + /// Validate minimum. + /// + [Test] + public void ValidateMinimum() + { + var n = new InverseGaussian(1.0, 2.0); + Assert.AreEqual(0.0, n.Minimum); + } + + /// + /// Validate maximum. + /// + [Test] + public void ValidateMaximum() + { + var n = new InverseGaussian(1.0, 2.0); + Assert.AreEqual(Double.PositiveInfinity, n.Maximum); + } + + /// + /// Validate entropy throws not supported exception + /// + /// Mu parameter. + /// Lambda parameter. + [TestCase(1.0, 1.0)] + public void ValidateEntropyFailsWithNotSupported(double mu, double lambda) + { + var n = new InverseGaussian(mu, lambda); + Assert.Throws(() => { var x = n.Entropy; }); + } + + /// + /// Validate skewness. + /// + /// Mu parameter. + /// Lambda parameter. + /// Expected value. + [TestCase(0.100000, 1.500000, 0.774596669241483)] + [TestCase(1.500000, 2.500000, 2.323790007724450)] + [TestCase(1.500000, 5.500000, 1.566698903601281)] + [TestCase(2.500000, 0.100000, 15.0)] + [TestCase(2.500000, 1.500000, 3.872983346207417)] + [TestCase(5.500000, 2.500000, 4.449719092257398)] + [TestCase(5.500000, 5.500000, 3.0)] + public void ValidateSkewness(double mu, double lambda, double skewness) + { + var n = new InverseGaussian(mu, lambda); + AssertHelpers.AlmostEqualRelative(skewness, n.Skewness, precision); + } + + /// + /// Validate density. + /// + /// Mu parameter. + /// Lambda parameter. + /// The location at which to compute the density. + /// Expected value. + [TestCase(1.500000, 2.500000, 0.500000, 0.587321148416470)] + [TestCase(1.500000, 2.500000, 0.800000, 0.627284170789435)] + [TestCase(2.500000, 2.500000, 0.500000, 0.360208446721537)] + [TestCase(2.500000, 2.500000, 0.800000, 0.428023216678204)] + public void ValidateDensity(double mu, double lambda, double x, double p) + { + var n = new InverseGaussian(mu, lambda); + AssertHelpers.AlmostEqualRelative(p, n.Density(x), precision); + AssertHelpers.AlmostEqualRelative(p, InverseGaussian.PDF(mu, lambda, x), precision); + } + + /// + /// Validate density log. + /// + /// Mu parameter. + /// Sigma value. + /// The location at which to compute the log density. + /// Expected value. + [TestCase(1.500000, 0.100000, 0.100000, 0.948091004233817)] + [TestCase(1.500000, 1.500000, 0.800000, -0.585657318845943)] + [TestCase(2.500000, 0.100000, 0.800000, -1.764415752730381)] + [TestCase(2.500000, 1.500000, 0.100000, -4.174328339659523)] + [TestCase(2.500000, 1.500000, 0.500000, -0.636485208310673)] + public void ValidateDensityLn(double mu, double lambda, double x, double p) + { + var n = new InverseGaussian(mu, lambda); + AssertHelpers.AlmostEqualRelative(p, n.DensityLn(x), precision); + AssertHelpers.AlmostEqualRelative(p, InverseGaussian.PDFLn(mu, lambda, x), precision); + } + + /// + /// Validate cumulative distribution. + /// + /// Mu parameter. + /// Sigma value. + /// The location at which to compute the density. + /// Expected value. + [TestCase(1.500000, 2.500000, 0.100000, 0.000002882116410867029)] + [TestCase(2.500000, 1.500000, 0.100000, 0.0001938001952257318)] + [TestCase(2.500000, 1.500000, 0.500000, 0.145457623130791)] + [TestCase(2.500000, 2.500000, 0.100000, 0.000001529605202470422)] + [TestCase(2.500000, 2.500000, 0.800000, 0.187168922781367)] + public void ValidateCumulativeDistribution(double mu, double lambda, double x, double f) + { + var n = new InverseGaussian(mu, lambda); + AssertHelpers.AlmostEqualRelative(f, n.CumulativeDistribution(x), precision); + AssertHelpers.AlmostEqualRelative(f, InverseGaussian.CDF(mu, lambda, x), precision); + } + + /// + /// Validate inverse cumulative distribution. + /// + /// Mu parameter. + /// Sigma value. + /// The location at which to compute the inverse cumulative distribution function. + /// Expected value. + [TestCase(1.500000, 1.500000, 0.100000, 0.356437063090717)] + [TestCase(1.500000, 1.500000, 0.500000, 1.013761958542859)] + [TestCase(2.500000, 2.500000, 0.500000, 1.689603264238098)] + [TestCase(2.500000, 2.500000, 0.800000, 3.619719792074686)] + public void ValidateInverseCumulativeDistribution(double mu, double lambda, double probability, double f) + { + var n = new InverseGaussian(mu, lambda); + AssertHelpers.AlmostEqualRelative(f, InverseGaussian.ICDF(mu, lambda, probability), precision); + AssertHelpers.AlmostEqualRelative(f, n.InvCDF(probability), precision); + } + + /// + /// Can sample static. + /// + [Test] + public void CanSampleStatic() + { + InverseGaussian.Sample(new Numerics.Random.MersenneTwister(100), 1.0, 1.0); + } + + /// + /// Can fill array with samples static. + /// + [Test] + public void CanFillSampleArrayStatic() + { + double[] samples = new double[100]; + InverseGaussian.Samples(new Numerics.Random.MersenneTwister(100), samples, 1.0, 1.0); + Assert.IsTrue(!samples.Any(x => x == 0)); + } + + /// + /// Can sample sequence static. + /// + [Test] + public void CanSampleSequenceStatic() + { + var ied = InverseGaussian.Samples(new Numerics.Random.MersenneTwister(100), 1.0, 1.0); + GC.KeepAlive(ied.Take(5).ToArray()); + } + + /// + /// Fail sample static with bad parameters. + /// + [Test] + public void FailSampleStatic() + { + Assert.That(() => { var d = InverseGaussian.Sample(new Numerics.Random.MersenneTwister(100), 0.0, -1.0); }, Throws.ArgumentException); + } + + /// + /// Fail filling array with samples with bad parameters static. + /// + [Test] + public void FailFillingSampleArrayStatic() + { + double[] samples = new double[100]; + Assert.That(() => { InverseGaussian.Samples(new Numerics.Random.MersenneTwister(100), samples, -1.0, 1.0); }, Throws.ArgumentException); + } + + /// + /// Fail sample sequence static with bad parameters. + /// + [Test] + public void FailSampleSequenceStatic() + { + Assert.That(() => { var ied = InverseGaussian.Samples(new Numerics.Random.MersenneTwister(100), 0.0, -1.0).First(); }, Throws.ArgumentException); + } + + /// + /// Can sample. + /// + [Test] + public void CanSample() + { + var n = new InverseGaussian(1.0, 2.0); + n.Sample(); + } + + /// + /// Can fill array with samples. + /// + [Test] + public void CanFillSampleArray() + { + double[] samples = new double[100]; + var n = new InverseGaussian(1.0, 2.0, new Numerics.Random.MersenneTwister(100)); + n.Samples(samples); + Assert.IsTrue(!samples.Any(x => x == 0)); + } + + /// + /// Can sample sequence. + /// + [Test] + public void CanSampleSequence() + { + var n = new InverseGaussian(1.0, 2.0); + var ied = n.Samples(); + GC.KeepAlive(ied.Take(5).ToArray()); + } + + /// + /// Can estimate distribution parameters. + /// + [TestCase(10.0, 0.1)] + [TestCase(1.5, 1.5)] + [TestCase(2.5, 2.5)] + [TestCase(2.5, 5.0)] + [TestCase(10.0, 50.0)] + public void CanEstimateParameters(double mu, double lambda) + { + var original = new InverseGaussian(mu, lambda, new Numerics.Random.MersenneTwister(100)); + var estimated = InverseGaussian.Estimate(original.Samples().Take(1000000)); + + AssertHelpers.AlmostEqualRelative(mu, estimated.Mu, 1); + AssertHelpers.AlmostEqualRelative(lambda, estimated.Lambda, 1); + } + } +} diff --git a/src/Numerics.Tests/DistributionTests/Continuous/TruncatedParetoTests.cs b/src/Numerics.Tests/DistributionTests/Continuous/TruncatedParetoTests.cs new file mode 100644 index 00000000..f5bae96a --- /dev/null +++ b/src/Numerics.Tests/DistributionTests/Continuous/TruncatedParetoTests.cs @@ -0,0 +1,446 @@ +// +// Math.NET Numerics, part of the Math.NET Project +// http://numerics.mathdotnet.com +// http://github.com/mathnet/mathnet-numerics +// +// Copyright (c) 2009-2019 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 +// restriction, including without limitation the rights to use, +// copy, modify, merge, publish, distribute, sublicense, and/or sell +// 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 +// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT +// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, +// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING +// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR +// OTHER DEALINGS IN THE SOFTWARE. +// + +using MathNet.Numerics.Distributions; +using NUnit.Framework; +using System; +using System.Linq; + +namespace MathNet.Numerics.UnitTests.DistributionTests.Continuous +{ + /// + /// TruncatedPareto distribution tests. + /// + [TestFixture, Category("Distributions")] + public class TruncatedParetoTests + { + private const int highPrecision = 12; + private const int lowPrecision = 6; + private const double highTruncation = 1E8; + + /// + /// Can create TruncatedPareto. + /// + /// Scale parameter. + /// Shape parameter. + /// /// Truncation parameter. + [TestCase(2.0, 2.0, 500.0)] + [TestCase(10.0, 2.0, 1000.0)] + [TestCase(100.0, 4.0, 10000.0)] + [TestCase(3.0, 10.0, 5.0)] + [TestCase(10.0, 100.0, 20.0)] + public void CanCreateTruncatedPareto(double scale, double shape, double truncation) + { + var n = new TruncatedPareto(scale, shape, truncation); + Assert.AreEqual(scale, n.Scale, highPrecision); + Assert.AreEqual(shape, n.Shape, highPrecision); + Assert.AreEqual(truncation, n.Truncation, highPrecision); + } + + /// + /// Can create TruncatedPareto with random source. + /// + [Test] + public void CanCreateTruncatedParetoWithRandomSource() + { + var randomSource = new Numerics.Random.MersenneTwister(100); + var n = new TruncatedPareto(10.0, 10.0, 1000, randomSource); + Assert.AreEqual(randomSource, n.RandomSource); + } + + /// + /// TruncatedPareto create fails with bad parameters. + /// + /// Scale parameter. + /// Shape parameter. + /// /// Truncation parameter. + [TestCase(1.0, -1.0, 15.0)] + [TestCase(1.0, 3.0, 0.5)] + [TestCase(-1.0, 2.0, 15.0)] + [TestCase(double.NaN, 1.0, 1.0)] + [TestCase(1.0, double.PositiveInfinity, 0.0)] + [TestCase(5.0, 2.0, double.PositiveInfinity)] + public void TruncatedParetoCreateFailsWithBadParameters(double scale, double shape, double truncation) + { + Assert.That(() => new TruncatedPareto(scale, shape, truncation), Throws.ArgumentException); + } + + /// + /// Validate IsValidParameterSet. + /// + /// Scale parameter. + /// Shape parameter. + /// /// Truncation parameter. + [TestCase(2.0, 2.0, 500, true)] + [TestCase(10.0, 10.0, 1000, true)] + [TestCase(5.0, 3.0, 1, false)] + public void ValidateIsValidParameterSet(double scale, double shape, double truncation, bool validity) + { + Assert.AreEqual(TruncatedPareto.IsValidParameterSet(scale, shape, truncation), validity); + } + + /// + /// Validate ToString. + /// + [Test] + public void ValidateToString() + { + var n = new TruncatedPareto(1d, 2d, 100d); + Assert.AreEqual("Truncated Pareto(Scale = 1, Shape = 2, Truncation = 100)", n.ToString()); + } + + /// + /// Validate mode throws not supported exception + /// + /// Scale value. + /// Shape parameter. + /// Truncation parameter. + [TestCase(10.0, 10.0, highTruncation)] + public void ValidateModeFailsWithNotSupported(double scale, double shape, double truncation) + { + var n = new TruncatedPareto(scale, shape, truncation); + Assert.Throws(() => { var x = n.Mode; }); + } + + /// + /// Validate median. + /// + /// Scale value. + /// Shape parameter. + /// Truncation parameter. + /// Expected value. + [TestCase(10.0, 10.0, highTruncation, 10.717734625362931)] + [TestCase(100, 3.5, 10000, 1.219013619375516e+02)] + [TestCase(1000, 5.5, 100000, 1.134312522193400e+03)] + public void ValidateMedian(double scale, double shape, double truncation, double median) + { + var n = new TruncatedPareto(scale, shape, truncation); + AssertHelpers.AlmostEqualRelative(median, n.Median, highPrecision); + } + + /// + /// Validate mean. + /// + /// Scale parameter. + /// Shape parameter. + /// Truncation parameter. + /// Expected value. + [TestCase(10.0, 10.0, highTruncation, 11.111111111111111)] + [TestCase(100, 3.5, 10000, 1.399986139998614e+02)] + public void ValidateMean(double scale, double shape, double truncation, double mean) + { + var n = new TruncatedPareto(scale, shape, truncation); + AssertHelpers.AlmostEqualRelative(mean, n.Mean, highPrecision); + } + + /// + /// Validate variance. + /// + /// Scale parameter. + /// Shape parameter. + /// Truncation parameter. + /// Expected value. + [TestCase(10.0, 10.0, highTruncation, 1.543209876543210)] + [TestCase(100, 3.5, 10000, 3.710390409118045e+03)] + [TestCase(1000, 5.5, 100000, 7.760125676537910e+04)] + public void ValidateVariance(double scale, double shape, double truncation, double variance) + { + var n = new TruncatedPareto(scale, shape, truncation); + AssertHelpers.AlmostEqualRelative(variance, n.Variance, highPrecision); + } + + /// + /// Validate standard deviation. + /// + /// Scale parameter. + /// Shape parameter. + /// Truncation parameter. + /// Expected value. + [TestCase(10.0, 10.0, highTruncation, 1.242259987499883)] + [TestCase(100, 3.5, 10000, 60.912974062329653)] + [TestCase(1000, 5.5, 100000, 2.785700212969427e+02)] + public void ValidateStandardDeviation(double scale, double shape, double truncation, double std) + { + var n = new TruncatedPareto(scale, shape, truncation); + AssertHelpers.AlmostEqualRelative(std, n.StdDev, highPrecision); + } + + /// + /// Validate minimum. + /// + [Test] + public void ValidateMinimum() + { + var n = new TruncatedPareto(1.0, 2.0, 500.0); + Assert.AreEqual(1.0, n.Minimum); + } + + /// + /// Validate maximum. + /// + [Test] + public void ValidateMaximum() + { + var n = new TruncatedPareto(1.0, 2.0, 500.0); + Assert.AreEqual(500.0, n.Maximum); + } + + /// + /// Validate entropy throws not supported exception + /// + /// Mu parameter. + /// Lambda parameter. + [TestCase(100, 3.5, 10000)] + public void ValidateEntropyFailsWithNotSupported(double scale, double shape, double truncation) + { + var n = new TruncatedPareto(scale, shape, truncation); + Assert.Throws(() => { var x = n.Entropy; }); + } + + /// + /// Validate GetMoments. + /// + /// Scale parameter. + /// Shape parameter. + /// Truncation parameter. + /// order of the moment. + /// Expected value. + [TestCase(10.0, 10.0, highTruncation, 1, 11.111111111111111)] + [TestCase(1.0, 6.0, 10000, 1, 1.2)] + [TestCase(100, 3.5, 10000, 1, 1.399986139998614e+02)] + [TestCase(100, 3.5, 10000, 2, 2.331000233100023e+04)] + [TestCase(100, 3.5, 10000, 3, 6.300000630000063e+06)] + [TestCase(1000, 5.5, 100000, 1, 1.222222221012222e+03)] + [TestCase(1000, 5.5, 100000, 2, 1.571428414301428e+06)] + [TestCase(1000, 5.5, 100000, 4, 3.663000000036630e+12)] + [TestCase(1000, 2.0, 100000, 2, 9.211261498125996e+06)] + [TestCase(1000, 3.0, 100000, 3, 1.381552437348865e+10)] + public void ValidateMoments(double scale, double shape, double truncation, int order, double value) + { + var n = new TruncatedPareto(scale, shape, truncation); + AssertHelpers.AlmostEqualRelative(value, n.GetMoment(order), highPrecision); + } + + /// + /// Validate skewness. + /// + /// Scale parameter. + /// Shape parameter. + /// Truncation parameter. + /// Expected value. + [TestCase(10.0, 10.0, highTruncation, 2.8110568859997356)] + public void ValidateSkewness(double scale, double shape, double truncation, double skewness) + { + var n = new TruncatedPareto(scale, shape, truncation); + AssertHelpers.AlmostEqualRelative(skewness, n.Skewness, highPrecision); + } + + /// + /// Validate density. + /// + /// Scale parameter. + /// Shape parameter. + /// Truncation parameter. + /// The location at which to compute the density. + /// Expected value. + [TestCase(1, 1, highTruncation, 1, 1)] + [TestCase(1, 1, highTruncation, 1.5, 4 / 9.0)] + [TestCase(1, 1, highTruncation, 5, 1 / 25.0)] + [TestCase(1, 1, highTruncation, 50, 1 / 2500.0)] + [TestCase(1, 4, highTruncation, 1, 4)] + [TestCase(1, 4, highTruncation, 1.5, 128 / 243.0)] + [TestCase(1, 4, highTruncation, 50, 1 / 78125000.0)] + [TestCase(3, 2, highTruncation, 3, 2 / 3.0)] + [TestCase(3, 2, highTruncation, 5, 18 / 125.0)] + [TestCase(25, 100, highTruncation, 50, 1.5777218104420236e-30)] + [TestCase(100, 25, highTruncation, 150, 6.6003546737276816e-6)] + [TestCase(100, 3.5, 10000, 1000, 1.106797291738662e-06)] + [TestCase(100, 3.5, 10000, 2000, 4.891399189920708e-08)] + [TestCase(1000, 5.5, 100000, 2000, 6.076698900882660e-05)] + public void ValidateDensity(double scale, double shape, double truncation, double x, double p) + { + var n = new TruncatedPareto(scale, shape, truncation); + AssertHelpers.AlmostEqualRelative(p, n.Density(x), lowPrecision); + AssertHelpers.AlmostEqualRelative(p, TruncatedPareto.PDF(scale, shape, truncation, x), lowPrecision); + } + + /// + /// Validate density log. + /// + /// Scale parameter. + /// Shape parameter. + /// Truncation parameter. + /// The location at which to compute the log density. + /// Expected value. + [TestCase(3, 2, highTruncation, 3, -0.405465108108164)] + [TestCase(100, 3.5, 10000, 1000, -13.714040035965924)] + [TestCase(100, 3.5, 10000, 2000, -16.833202348485678)] + public void ValidateDensityLn(double scale, double shape, double truncation, double x, double p) + { + var n = new TruncatedPareto(scale, shape, truncation); + AssertHelpers.AlmostEqualRelative(p, n.DensityLn(x), lowPrecision); + AssertHelpers.AlmostEqualRelative(p, TruncatedPareto.PDFLn(scale, shape, truncation, x), lowPrecision); + } + + /// + /// Validate cumulative distribution. + /// + /// Scale parameter. + /// Shape parameter. + /// Truncation parameter. + /// The location at which to compute the cumulative distribution function. + /// Expected value. + [TestCase(1.0, 1.0, highTruncation, 1.0, 0)] + [TestCase(7.0, 7.0, highTruncation, 10.0, 0.9176457)] + [TestCase(10.0, 10.0, highTruncation, 12.0, 0.83849441711015427)] + [TestCase(100, 3.5, 10000, 102, 0.066961862452387)] + [TestCase(100, 3.5, 10000, 1000, 0.999683872202370)] + [TestCase(100, 3.5, 10000, 2000, 0.999972149147496)] + [TestCase(1000, 5.5, 100000, 2000, 0.977902913097699)] + public void ValidateCumulativeDistribution(double scale, double shape, double truncation, double x, double expected) + { + var n = new TruncatedPareto(scale, shape, truncation); + AssertHelpers.AlmostEqualRelative(expected, n.CumulativeDistribution(x), highPrecision); + AssertHelpers.AlmostEqualRelative(expected, TruncatedPareto.CDF(scale, shape, truncation, x), highPrecision); + } + + /// + /// Validate inverse cumulative distribution. + /// + /// Scale parameter. + /// Shape parameter. + /// Truncation parameter. + /// The location at which to compute the inverse cumulative distribution function. + /// Expected value. + [TestCase(1.0, 1.0, highTruncation, 0, 1.0)] + [TestCase(7.0, 7.0, highTruncation, 0.9176457, 10.0)] + [TestCase(10.0, 10.0, highTruncation, 0.83849441711015427, 12.0)] + [TestCase(100, 3.5, 10000, 0.066961862452387, 102)] + [TestCase(100, 3.5, 10000, 0.999683872202370, 1000)] + [TestCase(100, 3.5, 10000, 0.999972149147496, 2000)] + [TestCase(1000, 5.5, 100000, 0.977902913097699, 2000)] + public void ValidateInverseCumulativeDistribution(double scale, double shape, double truncation, double p, double expected) + { + var n = new TruncatedPareto(scale, shape, truncation); + AssertHelpers.AlmostEqualRelative(expected, n.InvCDF(p), lowPrecision); + AssertHelpers.AlmostEqualRelative(expected, TruncatedPareto.ICDF(scale, shape, truncation, p), lowPrecision); + } + + /// + /// Can sample static. + /// + [Test] + public void CanSampleStatic() + { + TruncatedPareto.Sample(new Numerics.Random.MersenneTwister(100), 10.0, 10.0, 1000.0); + } + + /// + /// Can fill array with samples static. + /// + [Test] + public void CanFillSampleArrayStatic() + { + double[] samples = new double[100]; + TruncatedPareto.Samples(new Numerics.Random.MersenneTwister(100), samples, 10.0, 10.0, 1000.0); + Assert.IsTrue(!samples.Any(x => x == 0)); + } + + /// + /// Can sample sequence static. + /// + [Test] + public void CanSampleSequenceStatic() + { + var ied = TruncatedPareto.Samples(new Numerics.Random.MersenneTwister(100), 10.0, 10.0, 1000.0); + GC.KeepAlive(ied.Take(5).ToArray()); + } + + /// + /// Fail sample static with bad parameters. + /// + [Test] + public void FailSampleStatic() + { + Assert.That(() => { var d = TruncatedPareto.Sample(new Numerics.Random.MersenneTwister(100), 10.0, 10.0, 5.0); }, Throws.ArgumentException); + } + + /// + /// Fail filling array with samples with bad parameters static. + /// + [Test] + public void FailFillingSampleArrayStatic() + { + double[] samples = new double[100]; + Assert.That(() => { TruncatedPareto.Samples(new Numerics.Random.MersenneTwister(100), samples, 10.0, 10.0, 5.0); }, Throws.ArgumentException); + } + + /// + /// Fail sample sequence static with bad parameters. + /// + [Test] + public void FailSampleSequenceStatic() + { + Assert.That(() => { var ied = TruncatedPareto.Samples(new Numerics.Random.MersenneTwister(100), 10.0, 10.0, 5.0).First(); }, Throws.ArgumentException); + } + + /// + /// Can sample. + /// + [Test] + public void CanSample() + { + var n = new TruncatedPareto(10.0, 10.0, 1000.0); + n.Sample(); + } + + /// + /// Can fill array with samples. + /// + [Test] + public void CanFillSampleArray() + { + double[] samples = new double[100]; + var n = new TruncatedPareto(10.0, 10.0, 1000.0); + n.Samples(samples); + Assert.IsTrue(!samples.Any(x => x == 0)); + } + + /// + /// Can sample sequence. + /// + [Test] + public void CanSampleSequence() + { + var n = new TruncatedPareto(10.0, 10.0, 1000.0); + var ied = n.Samples(); + GC.KeepAlive(ied.Take(5).ToArray()); + } + } +} diff --git a/src/Numerics/Distributions/Burr.cs b/src/Numerics/Distributions/Burr.cs new file mode 100644 index 00000000..02162020 --- /dev/null +++ b/src/Numerics/Distributions/Burr.cs @@ -0,0 +1,429 @@ +// +// Math.NET Numerics, part of the Math.NET Project +// http://numerics.mathdotnet.com +// http://github.com/mathnet/mathnet-numerics +// +// Copyright (c) 2009-2019 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 +// restriction, including without limitation the rights to use, +// copy, modify, merge, publish, distribute, sublicense, and/or sell +// 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 +// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT +// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, +// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING +// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR +// OTHER DEALINGS IN THE SOFTWARE. +// + +using MathNet.Numerics.Properties; +using MathNet.Numerics.Random; +using System; +using System.Collections.Generic; + +namespace MathNet.Numerics.Distributions +{ + public class Burr : IContinuousDistribution + { + private System.Random _random; + + /// + /// Gets the scale (a) of the distribution. Range: a > 0. + /// + public double a { get; } + + /// + /// Gets the first shape parameter (c) of the distribution. Range: c > 0. + /// + public double c { get; } + + /// + /// Gets the second shape parameter (k) of the distribution. Range: k > 0. + /// + public double k { get; } + + /// + /// Initializes a new instance of the Burr Type XII class. + /// + /// The scale parameter a of the Burr distribution. Range: a > 0. + /// The first shape parameter c of the Burr distribution. Range: c > 0. + /// The second shape parameter k of the Burr distribution. Range: k > 0. + /// The random number generator which is used to draw random samples. Optional, can be null. + public Burr(double a, double c, double k, System.Random randomSource = null) + { + if (!IsValidParameterSet(a, c, k)) + { + throw new ArgumentException(Resources.InvalidDistributionParameters); + } + _random = randomSource ?? SystemRandomSource.Default; + this.a = a; + this.c = c; + this.k = k; + } + + /// + /// A string representation of the distribution. + /// + /// a string representation of the distribution. + public override string ToString() + { + return "Burr(a = " + a + ", c = " + c + ", k = " + k + ")"; + } + + /// + /// Tests whether the provided values are valid parameters for this distribution. + /// + /// The scale parameter a of the Burr distribution. Range: a > 0. + /// The first shape parameter c of the Burr distribution. Range: c > 0. + /// The second shape parameter k of the Burr distribution. Range: k > 0. + public static bool IsValidParameterSet(double a, double c, double k) + { + var allFinite = a.IsFinite() && c.IsFinite() && k.IsFinite(); + return allFinite && a > 0.0 && c > 0.0 && k > 0.0; + } + + /// + /// Gets the random number generator which is used to draw random samples. + /// + public System.Random RandomSource + { + get { return _random; } + set { _random = value ?? SystemRandomSource.Default; } + } + + /// + /// Gets the mean of the Burr distribution. + /// + public double Mean + { + get + { + return (1 / SpecialFunctions.Gamma(k)) * a * SpecialFunctions.Gamma(1 + 1 / c) * SpecialFunctions.Gamma(k - 1 / c); + } + } + + /// + /// Gets the variance of the Burr distribution. + /// + public double Variance + { + get + { + return (1 / SpecialFunctions.Gamma(k)) * Math.Pow(a, 2) * SpecialFunctions.Gamma(1 + 2 / c) * SpecialFunctions.Gamma(k - 2 / c) + - Math.Pow((1 / SpecialFunctions.Gamma(k)) * a * SpecialFunctions.Gamma(1 + 1 / c) * SpecialFunctions.Gamma(k - 1 / c), 2); + } + } + + /// + /// Gets the standard deviation of the Burr distribution. + /// + public double StdDev + { + get + { + return Math.Sqrt(Variance); + } + } + + /// + /// Gets the mode of the Burr distribution. + /// + public double Mode + { + get + { + return a * Math.Pow((c - 1) / (c * k + 1), 1 / c); + } + } + + /// + /// Gets the minimum of the Burr distribution. + /// + public double Minimum + { + get { return 0.0; } + } + + /// + /// Gets the maximum of the Burr distribution. + /// + public double Maximum + { + get { return double.PositiveInfinity; } + } + + /// + /// Gets the entropy of the Burr distribution (currently not supported). + /// + public double Entropy + { + get { throw new NotSupportedException(); } + } + + /// + /// Gets the skewness of the Burr distribution. + /// + public double Skewness + { + get + { + var mean = Mean; + var variance = Variance; + var std = StdDev; + return (GetMoment(3) - 3 * mean * variance - mean * mean * mean) / (std * std * std); + } + } + + /// + /// Gets the median of the Burr distribution. + /// + public double Median + { + get + { + return a * Math.Pow(Math.Pow(2, 1 / k) - 1, 1 / c); + } + } + + /// + /// Generates a sample from the Burr distribution. + /// + /// a sample from the distribution. + public double Sample() + { + return SampleUnchecked(_random, a, c, k); + } + + /// + /// Fills an array with samples generated from the distribution. + /// + /// The array to fill with the samples. + public void Samples(double[] values) + { + SamplesUnchecked(_random, values, a, c, k); + } + + /// + /// Generates a sequence of samples from the Burr distribution. + /// + /// a sequence of samples from the distribution. + public IEnumerable Samples() + { + return SamplesUnchecked(_random, a, c, k); + } + + /// + /// Generates a sample from the Burr distribution. + /// + /// The random number generator to use. + /// The scale parameter a of the Burr distribution. Range: a > 0. + /// The first shape parameter c of the Burr distribution. Range: c > 0. + /// The second shape parameter k of the Burr distribution. Range: k > 0. + /// a sample from the distribution. + public static double Sample(System.Random rnd, double a, double c, double k) + { + if (!IsValidParameterSet(a, c, k)) + { + throw new ArgumentException(Resources.InvalidDistributionParameters); + } + return SampleUnchecked(rnd, a, c, k); + } + + /// + /// Fills an array with samples generated from the distribution. + /// + /// The random number generator to use. + /// The array to fill with the samples. + /// The scale parameter a of the Burr distribution. Range: a > 0. + /// The first shape parameter c of the Burr distribution. Range: c > 0. + /// The second shape parameter k of the Burr distribution. Range: k > 0. + public static void Samples(System.Random rnd, double[] values, double a, double c, double k) + { + if (!IsValidParameterSet(a, c, k)) + { + throw new ArgumentException(Resources.InvalidDistributionParameters); + } + SamplesUnchecked(rnd, values, a, c, k); + } + + /// + /// Generates a sequence of samples from the Burr distribution. + /// + /// The random number generator to use. + /// The scale parameter a of the Burr distribution. Range: a > 0. + /// The first shape parameter c of the Burr distribution. Range: c > 0. + /// The second shape parameter k of the Burr distribution. Range: k > 0. + /// a sequence of samples from the distribution. + public static IEnumerable Samples(System.Random rnd, double a, double c, double k) + { + if (!IsValidParameterSet(a, c, k)) + { + throw new ArgumentException(Resources.InvalidDistributionParameters); + } + return SamplesUnchecked(rnd, a, c, k); + } + + internal static double SampleUnchecked(System.Random rnd, double a, double c, double k) + { + var k_inv = 1 / k; + var c_inv = 1 / c; + double u = rnd.NextDouble(); + return a * Math.Pow(Math.Pow(1 - u, -k_inv) - 1, c_inv); + } + + internal static void SamplesUnchecked(System.Random rnd, double[] values, double a, double c, double k) + { + if (values.Length == 0) + { + return; + } + var k_inv = 1 / k; + var c_inv = 1 / c; + double[] u = rnd.NextDoubles(values.Length); + + for (var j = 0; j < values.Length; ++j) + { + values[j] = a * Math.Pow(Math.Pow(1 - u[j], -k_inv) - 1, c_inv); + } + } + + internal static IEnumerable SamplesUnchecked(System.Random rnd, double a, double c, double k) + { + while (true) + { + yield return SampleUnchecked(rnd, a, c, k); + } + } + + /// + /// Gets the n-th raw moment of the distribution. + /// + /// The order (n) of the moment. Range: n ≥ 1. + /// the n-th moment of the distribution. + public double GetMoment(double n) + { + if (n > k * c) + { + throw new ArgumentException(Resources.ArgumentParameterSetInvalid); + } + var lambda_n = (n / c) * SpecialFunctions.Gamma(n / c) * SpecialFunctions.Gamma(k - n / c); + return Math.Pow(a, n) * lambda_n / SpecialFunctions.Gamma(k); + } + + /// + /// Computes the probability density of the distribution (PDF) at x, i.e. ∂P(X ≤ x)/∂x. + /// + /// The location at which to compute the density. + /// the density at . + /// + public double Density(double x) + { + return DensityImpl(a, c, k, x); + } + + /// + /// Computes the log probability density of the distribution (lnPDF) at x, i.e. ln(∂P(X ≤ x)/∂x). + /// + /// The location at which to compute the log density. + /// the log density at . + /// + public double DensityLn(double x) + { + return DensityLnImpl(a, c, k, 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 cumulative distribution at location . + /// + public double CumulativeDistribution(double x) + { + return CumulativeDistributionImpl(a, c, k, x); + } + + /// + /// Computes the probability density of the distribution (PDF) at x, i.e. ∂P(X ≤ x)/∂x. + /// + /// The scale parameter a of the Burr distribution. Range: a > 0. + /// The first shape parameter c of the Burr distribution. Range: c > 0. + /// The second shape parameter k of the Burr distribution. Range: k > 0. + /// The location at which to compute the density. + /// the density at . + /// + public static double PDF(double a, double c, double k, double x) + { + if (!IsValidParameterSet(a, c, k)) + { + throw new ArgumentException(Resources.InvalidDistributionParameters); + } + return DensityImpl(a, c, k, x); + } + + /// + /// Computes the log probability density of the distribution (lnPDF) at x, i.e. ln(∂P(X ≤ x)/∂x). + /// + /// The scale parameter a of the Burr distribution. Range: a > 0. + /// The first shape parameter c of the Burr distribution. Range: c > 0. + /// The second shape parameter k of the Burr distribution. Range: k > 0. + /// The location at which to compute the log density. + /// the log density at . + /// + public static double PDFLn(double a, double c, double k, double x) + { + if (!IsValidParameterSet(a, c, k)) + { + throw new ArgumentException(Resources.InvalidDistributionParameters); + } + return DensityLnImpl(a, c, k, x); + } + + /// + /// Computes the cumulative distribution (CDF) of the distribution at x, i.e. P(X ≤ x). + /// + /// The scale parameter a of the Burr distribution. Range: a > 0. + /// The first shape parameter c of the Burr distribution. Range: c > 0. + /// The second shape parameter k of the Burr distribution. Range: k > 0. + /// The location at which to compute the cumulative distribution function. + /// the cumulative distribution at location . + /// + public static double CDF(double a, double c, double k, double x) + { + if (!IsValidParameterSet(a, c, k)) + { + throw new ArgumentException(Resources.InvalidDistributionParameters); + } + return CumulativeDistributionImpl(a, c, k, x); + } + + internal static double DensityImpl(double a, double c, double k, double x) + { + var numerator = (k * c / a) * Math.Pow(x / a, c - 1); + var denominator = Math.Pow(1 + Math.Pow(x / a, c), k + 1); + return numerator / denominator; + } + + internal static double DensityLnImpl(double a, double c, double k, double x) + { + return Math.Log(DensityImpl(a, c, k, x)); + } + + internal static double CumulativeDistributionImpl(double a, double c, double k, double x) + { + var denominator = Math.Pow(1 + Math.Pow(x / a, c), k); + return 1 - 1 / denominator; + } + } +} diff --git a/src/Numerics/Distributions/InverseGaussian.cs b/src/Numerics/Distributions/InverseGaussian.cs new file mode 100644 index 00000000..4355ccac --- /dev/null +++ b/src/Numerics/Distributions/InverseGaussian.cs @@ -0,0 +1,449 @@ +// +// Math.NET Numerics, part of the Math.NET Project +// http://numerics.mathdotnet.com +// http://github.com/mathnet/mathnet-numerics +// +// Copyright (c) 2009-2019 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 +// restriction, including without limitation the rights to use, +// copy, modify, merge, publish, distribute, sublicense, and/or sell +// 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 +// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT +// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, +// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING +// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR +// OTHER DEALINGS IN THE SOFTWARE. +// + +using MathNet.Numerics.Properties; +using MathNet.Numerics.Random; +using MathNet.Numerics.Statistics; +using System; +using System.Collections.Generic; +using System.Linq; + +namespace MathNet.Numerics.Distributions +{ + public class InverseGaussian : IContinuousDistribution + { + private System.Random _random; + + /// + /// Gets the mean (μ) of the distribution. Range: μ > 0. + /// + public double Mu { get; } + + /// + /// Gets the shape (λ) of the distribution. Range: λ > 0. + /// + public double Lambda { get; } + + /// + /// Initializes a new instance of the InverseGaussian class. + /// + /// The mean (μ) of the distribution. Range: μ > 0. + /// The shape (λ) of the distribution. Range: λ > 0. + /// The random number generator which is used to draw random samples. + public InverseGaussian(double mu, double lambda, System.Random randomSource = null) + { + if (!IsValidParameterSet(mu, lambda)) + { + throw new ArgumentException(Resources.InvalidDistributionParameters); + } + _random = randomSource ?? SystemRandomSource.Default; + Mu = mu; + Lambda = lambda; + } + + /// + /// A string representation of the distribution. + /// + /// a string representation of the distribution. + public override string ToString() + { + return "InverseGaussian(μ = " + Mu + ", λ = " + Lambda + ")"; + } + + /// + /// Tests whether the provided values are valid parameters for this distribution. + /// + /// The mean (μ) of the distribution. Range: μ > 0. + /// The shape (λ) of the distribution. Range: λ > 0. + public static bool IsValidParameterSet(double mu, double lambda) + { + var allFinite = mu.IsFinite() && lambda.IsFinite(); + return allFinite && mu > 0.0 && lambda > 0.0; + } + + /// + /// Gets the random number generator which is used to draw random samples. + /// + public System.Random RandomSource + { + get { return _random; } + set { _random = value ?? SystemRandomSource.Default; } + } + + /// + /// Gets the mean of the Inverse Gaussian distribution. + /// + public double Mean + { + get + { + return Mu; + } + } + + /// + /// Gets the variance of the Inverse Gaussian distribution. + /// + public double Variance + { + get + { + return Math.Pow(Mu, 3) / Lambda; + } + } + + /// + /// Gets the standard deviation of the Inverse Gaussian distribution. + /// + public double StdDev + { + get + { + return Math.Sqrt(Variance); + } + } + + /// + /// Gets the median of the Inverse Gaussian distribution. + /// No closed form analytical expression exists, so this value is approximated numerically and can throw an exception. + /// + public double Median + { + get { return InvCDF(0.5); } + } + + /// + /// Gets the minimum of the Inverse Gaussian distribution. + /// + public double Minimum + { + get { return 0.0; } + } + + /// + /// Gets the maximum of the Inverse Gaussian distribution. + /// + public double Maximum + { + get { return double.PositiveInfinity; } + } + + /// + /// Gets the skewness of the Inverse Gaussian distribution. + /// + public double Skewness + { + get { return 3 * Math.Sqrt(Mu / Lambda); } + } + + /// + /// Gets the kurtosis of the Inverse Gaussian distribution. + /// + public double Kurtosis + { + get { return 15 * Mu / Lambda; } + } + + /// + /// Gets the mode of the Inverse Gaussian distribution. + /// + public double Mode + { + get { return Mu * (Math.Sqrt(1 + (9 * Mu * Mu) / (4 * Lambda * Lambda)) - (3 * Mu) / (2 * Lambda)); } + } + + /// + /// Gets the entropy of the Inverse Gaussian distribution (currently not supported). + /// + public double Entropy + { + get { throw new NotSupportedException(); } + } + + /// + /// Generates a sample from the inverse Gaussian distribution. + /// + /// a sample from the distribution. + public double Sample() + { + return SampleUnchecked(_random, Mu, Lambda); + } + + /// + /// Fills an array with samples generated from the distribution. + /// + /// The array to fill with the samples. + public void Samples(double[] values) + { + SamplesUnchecked(_random, values, Mu, Lambda); + } + + /// + /// Generates a sequence of samples from the inverse Gaussian distribution. + /// + /// a sequence of samples from the distribution. + public IEnumerable Samples() + { + return SamplesUnchecked(_random, Mu, Lambda); + } + + /// + /// Generates a sample from the inverse Gaussian distribution. + /// + /// The random number generator to use. + /// The mean (μ) of the distribution. Range: μ > 0. + /// The shape (λ) of the distribution. Range: λ > 0. + /// a sample from the distribution. + public static double Sample(System.Random rnd, double mu, double lambda) + { + if (!IsValidParameterSet(mu, lambda)) + { + throw new ArgumentException(Resources.InvalidDistributionParameters); + } + return SampleUnchecked(rnd, mu, lambda); + } + + /// + /// Fills an array with samples generated from the distribution. + /// + /// The random number generator to use. + /// The array to fill with the samples. + /// The mean (μ) of the distribution. Range: μ > 0. + /// The shape (λ) of the distribution. Range: λ > 0. + public static void Samples(System.Random rnd, double[] values, double mu, double lambda) + { + if (!IsValidParameterSet(mu, lambda)) + { + throw new ArgumentException(Resources.InvalidDistributionParameters); + } + SamplesUnchecked(rnd, values, mu, lambda); + } + + /// + /// Generates a sequence of samples from the Burr distribution. + /// + /// The random number generator to use. + /// The mean (μ) of the distribution. Range: μ > 0. + /// The shape (λ) of the distribution. Range: λ > 0. + /// a sequence of samples from the distribution. + public static IEnumerable Samples(System.Random rnd, double mu, double lambda) + { + if (!IsValidParameterSet(mu, lambda)) + { + throw new ArgumentException(Resources.InvalidDistributionParameters); + } + return SamplesUnchecked(rnd, mu, lambda); + } + + internal static double SampleUnchecked(System.Random rnd, double mu, double lambda) + { + double v = MathNet.Numerics.Distributions.Normal.Sample(rnd, 0, 1); + double test = rnd.NextDouble(); + return InverseGaussianSampleImpl(mu, lambda, v, test); + } + + internal static void SamplesUnchecked(System.Random rnd, double[] values, double mu, double lambda) + { + if (values.Length == 0) + { + return; + } + double[] v = new double[values.Length]; + MathNet.Numerics.Distributions.Normal.Samples(rnd, v, 0, 1); + double[] test = rnd.NextDoubles(values.Length); + for (var j = 0; j < values.Length; ++j) + { + values[j] = InverseGaussianSampleImpl(mu, lambda, v[j], test[j]); + } + } + + internal static IEnumerable SamplesUnchecked(System.Random rnd, double mu, double lambda) + { + while (true) + { + yield return SampleUnchecked(rnd, mu, lambda); + } + } + + internal static double InverseGaussianSampleImpl(double mu, double lambda, double normalSample, double uniformSample) + { + double y = normalSample * normalSample; + double x = mu + (mu * mu * y) / (2 * lambda) - (mu / (2 * lambda)) * Math.Sqrt(4 * mu * lambda * y + mu * mu * y * y); + if (uniformSample <= mu / (mu + x)) + return x; + else + return mu * mu / x; + } + + /// + /// Computes the probability density of the distribution (PDF) at x, i.e. ∂P(X ≤ x)/∂x. + /// + /// The location at which to compute the density. + /// the density at . + /// + public double Density(double x) + { + return DensityImpl(Mu, Lambda, x); + } + + /// + /// Computes the log probability density of the distribution (lnPDF) at x, i.e. ln(∂P(X ≤ x)/∂x). + /// + /// The location at which to compute the log density. + /// the log density at . + /// + public double DensityLn(double x) + { + return DensityLnImpl(Mu, Lambda, 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 cumulative distribution at location . + /// + public double CumulativeDistribution(double x) + { + return CumulativeDistributionImpl(Mu, Lambda, x); + } + + /// + /// Computes the inverse cumulative distribution (CDF) of the distribution at p, i.e. solving for P(X ≤ x) = p. + /// + /// The location at which to compute the inverse cumulative distribution function. + /// the inverse cumulative distribution at location . + public double InvCDF(double p) + { + Func equationToSolve = (x) => CumulativeDistribution(x) - p; + if (RootFinding.NewtonRaphson.TryFindRoot(equationToSolve, Density, Mode, 0, double.PositiveInfinity, 1e-8, 100, out double quantile)) + return quantile; + else + throw new NonConvergenceException(Resources.NumericalEstimationFailed); + } + + /// + /// Computes the probability density of the distribution (PDF) at x, i.e. ∂P(X ≤ x)/∂x. + /// + /// The mean (μ) of the distribution. Range: μ > 0. + /// The shape (λ) of the distribution. Range: λ > 0. + /// The location at which to compute the density. + /// the density at . + /// + public static double PDF(double mu, double lambda, double x) + { + if (!IsValidParameterSet(mu, lambda)) + { + throw new ArgumentException(Resources.InvalidDistributionParameters); + } + return DensityImpl(mu, lambda, x); + } + + /// + /// Computes the log probability density of the distribution (lnPDF) at x, i.e. ln(∂P(X ≤ x)/∂x). + /// + /// The mean (μ) of the distribution. Range: μ > 0. + /// The shape (λ) of the distribution. Range: λ > 0. + /// The location at which to compute the log density. + /// the log density at . + /// + public static double PDFLn(double mu, double lambda, double x) + { + if (!IsValidParameterSet(mu, lambda)) + { + throw new ArgumentException(Resources.InvalidDistributionParameters); + } + return DensityLnImpl(mu, lambda, x); + } + + /// + /// Computes the cumulative distribution (CDF) of the distribution at x, i.e. P(X ≤ x). + /// + /// The mean (μ) of the distribution. Range: μ > 0. + /// The shape (λ) of the distribution. Range: λ > 0. + /// The location at which to compute the cumulative distribution function. + /// the cumulative distribution at location . + /// + public static double CDF(double mu, double lambda, double x) + { + if (!IsValidParameterSet(mu, lambda)) + { + throw new ArgumentException(Resources.InvalidDistributionParameters); + } + return CumulativeDistributionImpl(mu, lambda, x); + } + + /// + /// Computes the inverse cumulative distribution (CDF) of the distribution at p, i.e. solving for P(X ≤ x) = p. + /// + /// The mean (μ) of the distribution. Range: μ > 0. + /// The shape (λ) of the distribution. Range: λ > 0. + /// The location at which to compute the inverse cumulative distribution function. + /// the inverse cumulative distribution at location . + /// + public static double ICDF(double mu, double lambda, double p) + { + if (!IsValidParameterSet(mu, lambda)) + { + throw new ArgumentException(Resources.InvalidDistributionParameters); + } + var igDist = new InverseGaussian(mu, lambda); + return igDist.InvCDF(p); + } + + /// + /// Estimates the Inverse Gaussian parameters from sample data with maximum-likelihood. + /// + /// The samples to estimate the distribution parameters from. + /// The random number generator which is used to draw random samples. Optional, can be null. + /// An Inverse Gaussian distribution. + public static InverseGaussian Estimate(IEnumerable samples, System.Random randomSource = null) + { + var sampleVec = samples.ToArray(); + var muHat = sampleVec.Mean(); + var lambdahat = 1 / (1 / samples.HarmonicMean() - 1 / muHat); + return new InverseGaussian(muHat, lambdahat, randomSource); + } + + internal static double DensityImpl(double mu, double lambda, double x) + { + return Math.Sqrt(lambda / (2 * Math.PI * Math.Pow(x, 3))) * Math.Exp(-((lambda * Math.Pow(x - mu, 2)) / (2 * mu * mu * x))); + } + + internal static double DensityLnImpl(double mu, double lambda, double x) + { + return Math.Log(Math.Sqrt(lambda / (2 * Math.PI * Math.Pow(x, 3)))) - ((lambda * Math.Pow(x - mu, 2)) / (2 * mu * mu * x)); + } + + internal static double CumulativeDistributionImpl(double mu, double lambda, double x) + { + return Normal.CDF(0, 1, Math.Sqrt(lambda / x) * (x / mu - 1)) + Math.Exp(2 * lambda / mu) * Normal.CDF(0, 1, -Math.Sqrt(lambda / x) * (x / mu + 1)); + } + } +} diff --git a/src/Numerics/Distributions/TruncatedPareto.cs b/src/Numerics/Distributions/TruncatedPareto.cs new file mode 100644 index 00000000..62190733 --- /dev/null +++ b/src/Numerics/Distributions/TruncatedPareto.cs @@ -0,0 +1,466 @@ +// +// Math.NET Numerics, part of the Math.NET Project +// http://numerics.mathdotnet.com +// http://github.com/mathnet/mathnet-numerics +// +// Copyright (c) 2009-2019 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 +// restriction, including without limitation the rights to use, +// copy, modify, merge, publish, distribute, sublicense, and/or sell +// 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 +// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT +// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, +// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING +// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR +// OTHER DEALINGS IN THE SOFTWARE. +// + +using MathNet.Numerics.Properties; +using MathNet.Numerics.Random; +using System; +using System.Collections.Generic; + +namespace MathNet.Numerics.Distributions +{ + public class TruncatedPareto : IContinuousDistribution + { + private System.Random _random; + + /// + /// Initializes a new instance of the TruncatedPareto class. + /// + /// The scale (xm) of the distribution. Range: xm > 0. + /// The shape (α) of the distribution. Range: α > 0. + /// The truncation (T) of the distribution. Range: T > xm. + /// The random number generator which is used to draw random samples. + /// If or are non-positive or if T ≤ xm. + public TruncatedPareto(double scale, double shape, double truncation, System.Random randomSource = null) + { + if (!IsValidParameterSet(scale, shape, truncation)) + { + throw new ArgumentException(Resources.InvalidDistributionParameters); + } + _random = randomSource ?? SystemRandomSource.Default; + Scale = scale; + Shape = shape; + Truncation = truncation; + } + + /// + /// A string representation of the distribution. + /// + /// a string representation of the distribution. + public override string ToString() + { + return "Truncated Pareto(Scale = " + Scale + ", Shape = " + Shape + ", Truncation = " + Truncation + ")"; + } + + /// + /// Tests whether the provided values are valid parameters for this distribution. + /// + /// The scale (xm) of the distribution. Range: xm > 0. + /// The shape (α) of the distribution. Range: α > 0. + /// The truncation (T) of the distribution. Range: T > xm. + public static bool IsValidParameterSet(double scale, double shape, double truncation) + { + var allFinite = scale.IsFinite() && shape.IsFinite() && truncation.IsFinite(); + return allFinite && scale > 0.0 && shape > 0.0 && truncation > scale; + } + + /// + /// Gets the random number generator which is used to draw random samples. + /// + public System.Random RandomSource + { + get { return _random; } + set { _random = value ?? SystemRandomSource.Default; } + } + + /// + /// Gets the scale (xm) of the distribution. Range: xm > 0. + /// + public double Scale { get; } + + /// + /// Gets the shape (α) of the distribution. Range: α > 0. + /// + public double Shape { get; } + + /// + /// Gets the truncation (T) of the distribution. Range: T > 0. + /// + public double Truncation { get; } + + /// + /// Gets the n-th raw moment of the distribution. + /// + /// The order (n) of the moment. Range: n ≥ 1. + /// the n-th moment of the distribution. + public double GetMoment(int n) + { + double moment; + if (Shape.AlmostEqual(n)) + { + moment = ((Shape * Math.Pow(Scale, n)) / (1 - Math.Pow(Scale / Truncation, Shape))) * Math.Log(Truncation / Scale); + } + else + { + moment = ((Shape * Math.Pow(Scale, n)) / (Shape - n)) * ((1 - Math.Pow((Scale / Truncation), (Shape - n))) / (1 - Math.Pow(Scale / Truncation, Shape))); + } + + return moment; + } + + /// + /// Gets the mean of the truncated Pareto distribution. + /// + public double Mean + { + get + { + return GetMoment(1); + } + } + + /// + /// Gets the variance of the truncated Pareto distribution. + /// + public double Variance + { + get + { + return GetMoment(2) - Math.Pow(GetMoment(1), 2); + } + } + + /// + /// Gets the standard deviation of the truncated Pareto distribution. + /// + public double StdDev + { + get + { + return Math.Sqrt(Variance); + } + } + + /// + /// Gets the mode of the truncated Pareto distribution (not supported). + /// + public double Mode + { + get { throw new NotSupportedException(); } + } + + /// + /// Gets the minimum of the truncated Pareto distribution. + /// + public double Minimum + { + get { return Scale; } + } + + /// + /// Gets the maximum of the truncated Pareto distribution. + /// + public double Maximum + { + get { return Truncation; } + } + + /// + /// Gets the entropy of the truncated Pareto distribution (not supported). + /// + public double Entropy + { + get { throw new NotSupportedException(); } + } + + /// + /// Gets the skewness of the truncated Pareto distribution. + /// + public double Skewness + { + get + { + var mean = Mean; + var variance = Variance; + var std = StdDev; + return (GetMoment(3) - 3.0 * mean * variance - mean * mean * mean) / (std * std * std); + } + } + + /// + /// Gets the median of the truncated Pareto distribution. + /// + public double Median + { + get + { + return Scale * Math.Pow(1.0 - (1.0 / 2.0) * (1.0 - Math.Pow(Scale / Truncation, Shape)), -(1.0 / Shape)); + } + } + + /// + /// Generates a sample from the truncated Pareto distribution. + /// + /// a sample from the distribution. + public double Sample() + { + return SampleUnchecked(_random, Scale, Shape, Truncation); + } + + /// + /// Fills an array with samples generated from the distribution. + /// + /// The array to fill with the samples. + public void Samples(double[] values) + { + SamplesUnchecked(_random, values, Scale, Shape, Truncation); + } + + /// + /// Generates a sequence of samples from the truncated Pareto distribution. + /// + /// a sequence of samples from the distribution. + public IEnumerable Samples() + { + return SamplesUnchecked(_random, Scale, Shape, Truncation); + } + + /// + /// Generates a sample from the truncated Pareto distribution. + /// + /// The random number generator to use. + /// The scale (xm) of the distribution. Range: xm > 0. + /// The shape (α) of the distribution. Range: α > 0. + /// The truncation (T) of the distribution. Range: T > xm. + /// a sample from the distribution. + public static double Sample(System.Random rnd, double scale, double shape, double truncation) + { + if (!IsValidParameterSet(scale, shape, truncation)) + { + throw new ArgumentException(Resources.InvalidDistributionParameters); + } + return SampleUnchecked(rnd, scale, shape, truncation); + } + + /// + /// Fills an array with samples generated from the distribution. + /// + /// The random number generator to use. + /// The array to fill with the samples. + /// The scale (xm) of the distribution. Range: xm > 0. + /// The shape (α) of the distribution. Range: α > 0. + /// The truncation (T) of the distribution. Range: T > xm. + public static void Samples(System.Random rnd, double[] values, double scale, double shape, double truncation) + { + if (!IsValidParameterSet(scale, shape, truncation)) + { + throw new ArgumentException(Resources.InvalidDistributionParameters); + } + SamplesUnchecked(rnd, values, scale, shape, truncation); + } + + /// + /// Generates a sequence of samples from the truncated Pareto distribution. + /// + /// The random number generator to use. + /// The scale (xm) of the distribution. Range: xm > 0. + /// The shape (α) of the distribution. Range: α > 0. + /// The truncation (T) of the distribution. Range: T > xm. + /// a sequence of samples from the distribution. + public static IEnumerable Samples(System.Random rnd, double scale, double shape, double truncation) + { + if (!IsValidParameterSet(scale, shape, truncation)) + { + throw new ArgumentException(Resources.InvalidDistributionParameters); + } + return SamplesUnchecked(rnd, scale, shape, truncation); + } + + internal static double SampleUnchecked(System.Random rnd, double scale, double shape, double truncation) + { + double uniform = rnd.NextDouble(); + return InvCDFUncheckedImpl(scale, shape, truncation, uniform); + } + + internal static void SamplesUnchecked(System.Random rnd, double[] values, double scale, double shape, double truncation) + { + if (values.Length == 0) + { + return; + } + double[] uniforms = rnd.NextDoubles(values.Length); + for (var j = 0; j < values.Length; ++j) + { + values[j] = InvCDFUncheckedImpl(scale, shape, truncation, uniforms[j]); + } + } + + internal static IEnumerable SamplesUnchecked(System.Random rnd, double scale, double shape, double truncation) + { + while (true) + { + yield return SampleUnchecked(rnd, scale, shape, truncation); + } + } + + /// + /// Computes the probability density of the distribution (PDF) at x, i.e. ∂P(X ≤ x)/∂x. + /// + /// The location at which to compute the density. + /// the density at . + /// + public double Density(double x) + { + return DensityImpl(Scale, Shape, Truncation, x); + } + + /// + /// Computes the log probability density of the distribution (lnPDF) at x, i.e. ln(∂P(X ≤ x)/∂x). + /// + /// The location at which to compute the log density. + /// the log density at . + /// + public double DensityLn(double x) + { + return DensityLnImpl(Scale, Shape, Truncation, 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 cumulative distribution at location . + /// + public double CumulativeDistribution(double x) + { + return CumulativeDistributionImpl(Scale, Shape, Truncation, x); + } + + /// + /// Computes the inverse cumulative distribution (CDF) of the distribution at p, i.e. solving for P(X ≤ x) = p. + /// + /// The location at which to compute the inverse cumulative distribution function. + /// the inverse cumulative distribution at location . + public double InvCDF(double p) + { + return InvCDFUncheckedImpl(Scale, Shape, Truncation, p); + } + + /// + /// Computes the inverse cumulative distribution (CDF) of the distribution at p, i.e. solving for P(X ≤ x) = p. + /// + /// The scale (xm) of the distribution. Range: xm > 0. + /// The shape (α) of the distribution. Range: α > 0. + /// The truncation (T) of the distribution. Range: T > xm. + /// The location at which to compute the inverse cumulative distribution function. + /// the inverse cumulative distribution at location . + /// + public static double ICDF(double scale, double shape, double truncation, double p) + { + if (!IsValidParameterSet(scale, shape, truncation)) + { + throw new ArgumentException(Resources.InvalidDistributionParameters); + } + return InvCDFUncheckedImpl(scale, shape, truncation, p); + } + + /// + /// 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 truncation (T) of the distribution. Range: T > xm. + /// The location at which to compute the density. + /// the density at . + /// + public static double PDF(double scale, double shape, double truncation, double x) + { + if (!IsValidParameterSet(scale, shape, truncation)) + { + throw new ArgumentException(Resources.InvalidDistributionParameters); + } + return DensityImpl(scale, shape, truncation, x); + } + + /// + /// 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 truncation (T) of the distribution. Range: T > xm. + /// The location at which to compute the log density. + /// the log density at . + /// + public static double PDFLn(double scale, double shape, double truncation, double x) + { + if (!IsValidParameterSet(scale, shape, truncation)) + { + throw new ArgumentException(Resources.InvalidDistributionParameters); + } + return DensityLnImpl(scale, shape, truncation, x); + } + + /// + /// Computes the cumulative distribution (CDF) of the distribution at x, i.e. P(X ≤ x). + /// + /// The scale (xm) of the distribution. Range: xm > 0. + /// The shape (α) of the distribution. Range: α > 0. + /// The truncation (T) of the distribution. Range: T > xm. + /// The location at which to compute the cumulative distribution function. + /// the cumulative distribution at location . + /// + public static double CDF(double scale, double shape, double truncation, double x) + { + if (!IsValidParameterSet(scale, shape, truncation)) + { + throw new ArgumentException(Resources.InvalidDistributionParameters); + } + return CumulativeDistributionImpl(scale, shape, truncation, x); + } + + internal static double DensityImpl(double scale, double shape, double truncation, double x) + { + if (x < scale || x > truncation) + return 0; + else + return (shape * Math.Pow(scale, shape) * Math.Pow(x, -shape - 1)) / (1 - Math.Pow(scale / truncation, shape)); + } + + internal static double DensityLnImpl(double scale, double shape, double truncation, double x) + { + return Math.Log(DensityImpl(scale, shape, truncation, x)); + } + + internal static double CumulativeDistributionImpl(double scale, double shape, double truncation, double x) + { + if (x <= scale) + return 0; + else if (x >= truncation) + return 1; + else + return (1 - Math.Pow(scale, shape) * Math.Pow(x, -shape)) / (1 - Math.Pow(scale / truncation, shape)); + } + + internal static double InvCDFUncheckedImpl(double scale, double shape, double truncation, double p) + { + var numerator = p * Math.Pow(truncation, shape) - p * Math.Pow(scale, shape) - Math.Pow(truncation, shape); + var denominator = Math.Pow(truncation, shape) * Math.Pow(scale, shape); + return Math.Pow(-numerator / denominator, -(1 / shape)); + } + } +} diff --git a/src/Numerics/Precision.Comparison.cs b/src/Numerics/Precision.Comparison.cs index 912d8ab6..c437494b 100644 --- a/src/Numerics/Precision.Comparison.cs +++ b/src/Numerics/Precision.Comparison.cs @@ -674,5 +674,14 @@ namespace MathNet.Numerics return CompareToNumbersBetween(a, b, maxNumbersBetween) < 0; } + + /// + /// Checks if a given double values is finite, i.e. neither NaN nor inifnity + /// + /// The value to be checked fo finitenes. + public static bool IsFinite(this double value) + { + return !double.IsNaN(value) && !double.IsInfinity(value); + } } } diff --git a/src/Numerics/Properties/Resources.Designer.cs b/src/Numerics/Properties/Resources.Designer.cs index 2cd21dec..6556a755 100644 --- a/src/Numerics/Properties/Resources.Designer.cs +++ b/src/Numerics/Properties/Resources.Designer.cs @@ -853,6 +853,17 @@ namespace MathNet.Numerics.Properties } } + /// + /// Looks up a localized string similar to: Numerical estimation of the statistic has failed. + /// + public static string NumericalEstimationFailed + { + get + { + return ResourceManager.GetString("NumericalEstimationFailed", resourceCulture); + } + } + /// /// Looks up a localized string similar to The two arguments can't be compared (maybe they are part of a partial ordering?). /// diff --git a/src/Numerics/Properties/Resources.resx b/src/Numerics/Properties/Resources.resx index 14692289..6f8c7f4f 100644 --- a/src/Numerics/Properties/Resources.resx +++ b/src/Numerics/Properties/Resources.resx @@ -448,4 +448,7 @@ The sample size must be larger than the given degrees of freedom. + + Numerical estimation of the statistic has failed. The used solver did not succeed in finding a root. + \ No newline at end of file