diff --git a/src/FSharp/Distributions.fs b/src/FSharp/Distributions.fs index 674e48cd..d6f6a8a9 100644 --- a/src/FSharp/Distributions.fs +++ b/src/FSharp/Distributions.fs @@ -150,6 +150,14 @@ module Sample = let studentT location scale freedom (rng:System.Random) = StudentT.Sample(rng, location, scale, freedom) let studentTSeq location scale freedom (rng:System.Random) = StudentT.Samples(rng, location, scale, freedom) + /// Skew Generalized T with location (μ), scale (σ), skew (λ), kurtosis param (p) and kurtosis param (q). + let skewGeneralizedT location scale skew p q (rng:System.Random) = SkewedGeneralizedT.Sample(rng, location, scale, skew, p, q) + let skewGeneralizedTSeq location scale skew p q (rng:System.Random) = SkewedGeneralizedT.Samples(rng, location, scale, skew, p, q) + + /// Skew Generalized Error with location (μ), scale (σ), skew (λ) and kurtosis param (p). + let skewGeneralizedError location scale skew p (rng:System.Random) = SkewedGeneralizedError.Sample(rng, location, scale, skew, p) + let skewGeneralizedErrorSeq location scale skew p (rng:System.Random) = SkewedGeneralizedError.Samples(rng, location, scale, skew, p) + /// Weibull with shape (k) and scale (λ). let weibull shape scale (rng:System.Random) = Weibull.Sample(rng, shape, scale) let weibullSeq shape scale (rng:System.Random) = Weibull.Samples(rng, shape, scale) diff --git a/src/Numerics.Tests/DistributionTests/Continuous/SkewGeneralizedTTests.cs b/src/Numerics.Tests/DistributionTests/Continuous/SkewGeneralizedTTests.cs new file mode 100644 index 00000000..659e43da --- /dev/null +++ b/src/Numerics.Tests/DistributionTests/Continuous/SkewGeneralizedTTests.cs @@ -0,0 +1,261 @@ +// +// Math.NET Numerics, part of the Math.NET Project +// http://numerics.mathdotnet.com +// http://github.com/mathnet/mathnet-numerics +// +// Copyright (c) 2009-2016 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 System; +using System.Linq; +using MathNet.Numerics.Distributions; +using NUnit.Framework; + +namespace MathNet.Numerics.UnitTests.DistributionTests.Continuous +{ + /// + /// SkewedGeneralizedT distribution tests. + /// Reference values are from the R package sgt 2.0 (run on Microsoft R Open v3.5.2) + /// + [TestFixture, Category("Distributions")] + public class SkewedGeneralizedTTests + { + [Test] + public void CanCreateStandardSkewedGeneralizedT() + { + var n = new SkewedGeneralizedT(); + Assert.AreEqual(0.0, n.Location); + Assert.AreEqual(1.0, n.Scale); + Assert.AreEqual(0.0, n.Skew); + Assert.AreEqual(2.0, n.P); + Assert.AreEqual(double.PositiveInfinity, n.Q); + } + + [TestCase(0.0, 1.0, 0.0, 2.0, double.PositiveInfinity)] // Standard Normal distribution + [TestCase(5.0, 1.0, 0.0, 2.0, double.PositiveInfinity)] // Mean shifted Normal distribution + [TestCase(0.0, 2.0, 0.0, 2.0, double.PositiveInfinity)] // Scaled Normal distribution + [TestCase(0.0, 1.0, 0.9, 2.0, double.PositiveInfinity)] // Skewed Normal distribution + [TestCase(1.0, 1.5, 0.9, 2.0, double.PositiveInfinity)] // Mean shifted and scaled Skewed Normal distribution + [TestCase(0.0, 1.0, 0.9, 2.0, 1.1)] // Skewed Student T distribution + [TestCase(0.0, 1.0, 0.0, 2.0, 1.1)] // Student T distribution + [TestCase(0.0, 1.0, -0.3, 2.2, double.PositiveInfinity)] // Skewed Generalized Error Distribution + [TestCase(0.0, 1.0, 0.0, 2.2, double.PositiveInfinity)] // Generalized Error Distribution + [TestCase(0.0, 1.0, 0.0, double.PositiveInfinity, 1.0)] // Continuous Uniform + public void CanCreateSkewedGeneralizedT(double location, double scale, double skew, double p, double q) + { + var n = new SkewedGeneralizedT(location, scale, skew, p, q); + Assert.AreEqual(location, n.Location); + Assert.AreEqual(scale, n.Scale); + Assert.AreEqual(skew, n.Skew); + Assert.AreEqual(p, n.P); + Assert.AreEqual(q, n.Q); + } + + [TestCase(0.0, 1.0, 0.0, 2.0, 1.0)] // pq <= 2 + [TestCase(0.0, 1.0, 0.0, -2.0, -1.0)] // pq <= 2 and negative values + [TestCase(5.0, -1.0, 0.0, 2.0, double.PositiveInfinity)] // Negative scale + [TestCase(0.0, 2.0, 1.1, 2.0, double.PositiveInfinity)] // Invalid skew, too large + [TestCase(0.0, 1.0, -1.1, 2.0, double.PositiveInfinity)] // Invalid skew, too small + public void SkewedGeneralizedTCreateFailsWithBadParameters(double location, double scale, double skew, double p, double q) + { + Assert.That(() => new SkewedGeneralizedT(location, scale, skew, p, q), Throws.ArgumentException); + } + + [TestCase(0.0, 1.0, 0.0, 2.0, double.PositiveInfinity, 0.9123, 0.2631360242)] // Standard Normal distribution + [TestCase(5.0, 1.0, 0.0, 2.0, double.PositiveInfinity, 4.9123, 0.3974110362)] // Mean shifted Normal distribution + [TestCase(0.0, 2.0, 0.0, 2.0, double.PositiveInfinity, 0.9123, 0.1797618922)] // Scaled Normal distribution + [TestCase(0.0, 1.0, 0.9, 2.0, double.PositiveInfinity, 0.9123, 0.1958872375)] // Skewed Normal distribution + [TestCase(1.0, 1.5, 0.9, 2.0, double.PositiveInfinity, 0.9123, 0.2400040926)] // Mean shifted and scaled Skewed Normal distribution + [TestCase(0.0, 1.0, -0.9, 2.0, 5.0, 0.9123, 0.523647666)] // Skewed Student T distribution + [TestCase(0.0, 1.0, -0.9, 2.0, 5.0, -0.9123, 0.1799965988)] // Skewed Student T distribution + [TestCase(0.0, 1.0, 0.0, 2.0, 5.0, 0.9123, 0.2524160191)] // Student T distribution + [TestCase(0.0, 1.0, -0.3, 2.2, double.PositiveInfinity, 0.9123, 0.3323206895)] // Skewed Generalized Error Distribution + [TestCase(-1.5, 1.0, -0.3, 2.2, double.PositiveInfinity, 0.9123, 0.006297697854)] // Mean shifted Skewed Generalized Error Distribution + [TestCase(0.0, 1.0, 0.0, 2.2, double.PositiveInfinity, 0.9123, 0.2701962342)] // Generalized Error Distribution + [TestCase(0.0, 1.0, 0.0, double.PositiveInfinity, 1.0, 0.5123, 0.2886751346)] // Continuous Uniform + [TestCase(0.0, 1.0, 0.0, double.PositiveInfinity, 1.0, 0.6123, 0.2886751346)] // Continuous Uniform + [TestCase(0.0, 1.0, 0.77, 1.0, double.PositiveInfinity, 0.6123, 0.2016342715)] // Skewed Laplace + [TestCase(0.0, 1.0, 0.0, 1.0, double.PositiveInfinity, 0.6123, 0.2974537422)] // Laplace + [TestCase(0.9, 1.0, 0.0, 1.0, double.PositiveInfinity, 0.6123, 0.4707430703)] // Mean shifted Laplace + public void ValidateDensity(double location, double scale, double skew, double p, double q, double x, double d) + { + var n = new SkewedGeneralizedT(location, scale, skew, p, q); + var density = n.Density(x); + AssertHelpers.AlmostEqualRelative(d, density, 8); + } + + [TestCase(0.0, 1.0, 0.0, 2.0, double.PositiveInfinity, 0.9123, -1.335084178)] // Standard Normal distribution + [TestCase(5.0, 1.0, 0.0, 2.0, double.PositiveInfinity, 4.9123, -0.9227841782)] // Mean shifted Normal distribution + [TestCase(0.0, 2.0, 0.0, 2.0, double.PositiveInfinity, 0.9123, -1.716122125)] // Scaled Normal distribution + [TestCase(0.0, 1.0, 0.9, 2.0, double.PositiveInfinity, 0.9123, -1.630216104)] // Skewed Normal distribution + [TestCase(1.0, 1.5, 0.9, 2.0, double.PositiveInfinity, 0.9123, -1.427099303)] // Mean shifted and scaled Skewed Normal distribution + [TestCase(0.0, 1.0, -0.9, 2.0, 5.0, 0.9123, -0.646936214)] // Skewed Student T distribution + [TestCase(0.0, 1.0, -0.9, 2.0, 5.0, -0.9123, -1.714817324)] // Skewed Student T distribution + [TestCase(0.0, 1.0, 0.0, 2.0, 5.0, 0.9123, -1.376676683)] // Student T distribution + [TestCase(0.0, 1.0, -0.3, 2.2, double.PositiveInfinity, 0.9123, -1.101654844)] // Skewed Generalized Error Distribution + [TestCase(-1.5, 1.0, -0.3, 2.2, double.PositiveInfinity, 0.9123, -5.067571132)] // Mean shifted Skewed Generalized Error Distribution + [TestCase(0.0, 1.0, 0.0, 2.2, double.PositiveInfinity, 0.9123, -1.308606791)] // Generalized Error Distribution + [TestCase(0.0, 1.0, 0.0, double.PositiveInfinity, 1.0, -0.5123, -1.242453325)] // Continuous Uniform + [TestCase(0.0, 1.0, 0.0, double.PositiveInfinity, 1.0, -0.6123, -1.242453325)] // Continuous Uniform + [TestCase(0.0, 1.0, 0.77, 1.0, double.PositiveInfinity, 0.6123, -1.60129976)] // Skewed Laplace + [TestCase(0.0, 1.0, 0.0, 1.0, double.PositiveInfinity, 0.6123, -1.212496555)] // Laplace + [TestCase(0.9, 1.0, 0.0, 1.0, double.PositiveInfinity, 0.6123, -0.7534428322)] // Mean shifted Laplace + public void ValidateDensityLn(double location, double scale, double skew, double p, double q, double x, double d) + { + var n = new SkewedGeneralizedT(location, scale, skew, p, q); + var density = n.DensityLn(x); + AssertHelpers.AlmostEqualRelative(d, density, 8); + } + + [TestCase(0.0, 1.0, 0.0, 2.0, double.PositiveInfinity, 0.9123, 0.8191945928)] // Standard Normal distribution + [TestCase(5.0, 1.0, 0.0, 2.0, double.PositiveInfinity, 4.123, 0.190243319)] // Mean shifted Normal distribution + [TestCase(0.0, 2.0, 0.0, 2.0, double.PositiveInfinity, 0.9123, 0.6758589413)] // Scaled Normal distribution + [TestCase(0.0, 1.0, 0.9, 2.0, double.PositiveInfinity, 0.9123, 0.8216671619)] // Skewed Normal distribution + [TestCase(1.0, 1.5, 0.9, 2.0, double.PositiveInfinity, 0.9123, 0.5519973476)] // Mean shifted and scaled Skewed Normal distribution + [TestCase(0.0, 1.0, -0.9, 2.0, 5.0, 0.9123, 0.8297526431)] // Skewed Student T distribution + [TestCase(0.0, 1.0, -0.9, 2.0, 5.0, -0.9123, 0.1624618933)] // Skewed Student T distribution + [TestCase(0.0, 1.0, 0.0, 2.0, 5.0, 0.9123, 0.8341106883)] // Student T distribution + [TestCase(0.0, 1.0, -0.3, 2.2, double.PositiveInfinity, 0.9123, 0.8118701776)] // Skewed Generalized Error Distribution + [TestCase(-1.5, 1.0, -0.3, 2.2, double.PositiveInfinity, 0.9123, 0.9987893207)] // Mean shifted Skewed Generalized Error Distribution + [TestCase(0.0, 1.0, 0.0, 2.2, double.PositiveInfinity, 0.9123, 0.8140902875)] // Generalized Error Distribution + [TestCase(0.0, 1.0, 0.0, double.PositiveInfinity, 1.0, 0.5123, 0.6478882715)] // Continuous Uniform + [TestCase(0.0, 1.0, 0.0, double.PositiveInfinity, 1.0, 0.6123, 0.6767557849)] // Continuous Uniform + [TestCase(0.0, 1.0, 0.77, 1.0, double.PositiveInfinity, 0.6123, 0.8000467981)] // Skewed Laplace + [TestCase(0.0, 1.0, 0.0, 1.0, double.PositiveInfinity, 0.6123, 0.7896684418)] // Laplace + [TestCase(0.9, 1.0, 0.0, 1.0, double.PositiveInfinity, 0.6123, 0.3328656172)] // Mean shifted Laplace + public void ValidateCDF(double location, double scale, double skew, double p, double q, double x, double pr) + { + var n = new SkewedGeneralizedT(location, scale, skew, p, q); + var cpr = n.CumulativeDistribution(x); + AssertHelpers.AlmostEqualRelative(pr, cpr, 8); + } + + [TestCase(0.0, 1.0, 0.0, 2.0, double.PositiveInfinity, 0.9123, 1.355055108)] // Standard Normal distribution + [TestCase(5.0, 1.0, 0.0, 2.0, double.PositiveInfinity, 0.4123, 4.778367525)] // Mean shifted Normal distribution + [TestCase(0.0, 2.0, 0.0, 2.0, double.PositiveInfinity, 0.9123, 2.710110216)] // Scaled Normal distribution + [TestCase(0.0, 1.0, 0.9, 2.0, double.PositiveInfinity, 0.9123, 1.506912119)] // Skewed Normal distribution + [TestCase(1.0, 1.5, 0.9, 2.0, double.PositiveInfinity, 0.9123, 3.260368178)] // Mean shifted and scaled Skewed Normal distribution + [TestCase(0.0, 1.0, -0.9, 2.0, 5.0, 0.9123, 1.068716999)] // Skewed Student T distribution + [TestCase(0.0, 1.0, -0.9, 2.0, 5.0, 0.123, -1.159154003)] // Skewed Student T distribution + [TestCase(0.0, 1.0, 0.0, 2.0, 5.0, 0.9123, 1.304471125)] // Student T distribution + [TestCase(0.0, 1.0, -0.3, 2.2, double.PositiveInfinity, 0.9123, 1.275451275)] // Skewed Generalized Error Distribution + [TestCase(-1.5, 1.0, -0.3, 2.2, double.PositiveInfinity, 0.9123, -0.2245487247)] // Mean shifted Skewed Generalized Error Distribution + [TestCase(0.0, 1.0, 0.0, 2.2, double.PositiveInfinity, 0.9123, 1.363253899)] // Generalized Error Distribution + [TestCase(0.0, 1.0, 0.0, double.PositiveInfinity, 1.0, 0.5123, 0.04260844987)] // Continuous Uniform + [TestCase(0.0, 1.0, 0.0, double.PositiveInfinity, 1.0, 0.6123, 0.3890186114)] // Continuous Uniform + [TestCase(0.0, 1.0, 0.77, 1.0, double.PositiveInfinity, 0.6123, -0.04432801722)] // Skewed Laplace + [TestCase(0.0, 1.0, 0.0, 1.0, double.PositiveInfinity, 0.6123, 0.179871174)] // Laplace + [TestCase(0.9, 1.0, 0.0, 1.0, double.PositiveInfinity, 0.6123, 1.079871174)] // Mean shifted Laplace + public void ValidateInvCDF(double location, double scale, double skew, double p, double q, double quantile, double x) + { + var n = new SkewedGeneralizedT(location, scale, skew, p, q); + var xq = n.InverseCumulativeDistribution(quantile); + AssertHelpers.AlmostEqualRelative(x, xq, 8); + + AssertHelpers.AlmostEqualRelative(quantile, n.CumulativeDistribution(xq), 8); + } + + [TestCase(0, 1.0, 0.4123)] + [TestCase(1.5, 2.5, 0.5123)] + [TestCase(-0.5, 5, 0.6123)] + public void ValidateLaplaceDensityEquivalence(double location, double scale, double x) + { + var n = new SkewedGeneralizedT(location, scale, 0, 1, double.PositiveInfinity); + + var b = scale / Math.Sqrt(2.0); + var l = new Laplace(location, b); + + AssertHelpers.AlmostEqualRelative(l.Density(x), n.Density(x), 8); + AssertHelpers.AlmostEqualRelative(l.DensityLn(x), n.DensityLn(x), 8); + } + + [TestCase(0, 1.0, 0.4123)] + [TestCase(1.5, 2.5, 0.5123)] + [TestCase(-0.5, 5, 0.6123)] + public void ValidateNormalDensityEquivalence(double location, double scale, double x) + { + var sgt = new SkewedGeneralizedT(location, scale, 0, 2, double.PositiveInfinity); + var n = new Normal(location, scale); + + AssertHelpers.AlmostEqualRelative(n.Density(x), sgt.Density(x), 8); + AssertHelpers.AlmostEqualRelative(n.DensityLn(x), sgt.DensityLn(x), 8); + } + + [TestCase(0, 1, -0.1, 0.5123)] + [TestCase(0, 1, 0.1, 0.6123)] + public void ValidateSkewedNormalDistribution(double location, double scale, double skew, double x) + { + var sn = new SkewedGeneralizedT(location, scale, skew, 2, double.PositiveInfinity); + var n = new Normal(location, scale); + + var sp = sn.CumulativeDistribution(x); + var p = n.CumulativeDistribution(x); + + if (skew > 0) + Assert.IsTrue(sp > p); + else + Assert.IsTrue(sp < p); + } + + /// + /// Can sample static. + /// + [Test] + public void CanSampleStatic() + { + SkewedGeneralizedT.Sample(0.0, 1.0, 0.3, 2.2, 5.6); + } + + /// + /// Can sample sequence static. + /// + [Test] + public void CanSampleSequenceStatic() + { + var ied = SkewedGeneralizedT.Samples(0.0, 1.0, 0.3, 2.2, 5.6); + GC.KeepAlive(ied.Take(5).ToArray()); + } + + /// + /// Can sample. + /// + [Test] + public void CanSample() + { + var n = new SkewedGeneralizedT(); + n.Sample(); + } + + /// + /// Can sample sequence. + /// + [Test] + public void CanSampleSequence() + { + var n = new SkewedGeneralizedT(); + var ied = n.Samples(); + GC.KeepAlive(ied.Take(5).ToArray()); + } + } +} diff --git a/src/Numerics/Distributions/SkewedGeneralizedError.cs b/src/Numerics/Distributions/SkewedGeneralizedError.cs new file mode 100644 index 00000000..78f57fd7 --- /dev/null +++ b/src/Numerics/Distributions/SkewedGeneralizedError.cs @@ -0,0 +1,372 @@ +// +// 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 +{ + /// + /// Continuous Univariate Skewed Generalized Error Distribution (SGED). + /// Implements the univariate SSkewed Generalized Error Distribution. For details about this + /// distribution, see + /// + /// Wikipedia - Generalized Error Distribution. + /// It includes Laplace, Normal and Student-t distributions. + /// This is the distribution with q=Inf. + /// + /// This implementation is based on the R package dsgt and corresponding viginette, see + /// https://cran.r-project.org/web/packages/sgt/vignettes/sgt.pdf. Compared to that + /// implementation, the options for mean adjustment and variance adjustment are always true. + /// The location (μ) is the mean of the distribution. + /// The scale (σ) squared is the variance of the distribution. + /// + /// The distribution will use the by + /// default. Users can get/set the random number generator by using the + /// property. + /// The statistics classes will check all the incoming parameters + /// whether they are in the allowed range. + public class SkewedGeneralizedError : IContinuousDistribution + { + private System.Random _random; + + /// + /// Initializes a new instance of the SkewedGeneralizedError class. This is a generalized error distribution + /// with location=0.0, scale=1.0, skew=0.0 and p=2.0 (a standard normal distribution). + /// + public SkewedGeneralizedError() + { + _random = SystemRandomSource.Default; + Location = 0.0; + Scale = 1.0; + Skew = 0.0; + P = 2.0; + } + + /// + /// Initializes a new instance of the SkewedGeneralizedT class with a particular location, scale, skew + /// and kurtosis parameters. Different parameterizations result in different distributions. + /// + /// The location (μ) of the distribution. + /// The scale (σ) of the distribution. Range: σ > 0. + /// The skew, 1 > λ > -1 + /// Parameter that controls kurtosis. Range: p > 0 + public SkewedGeneralizedError(double location, double scale, double skew, double p) + { + if (!IsValidParameterSet(location, scale, skew, p)) + { + throw new ArgumentException(Resources.InvalidDistributionParameters); + } + + _random = SystemRandomSource.Default; + Location = location; + Scale = scale; + Skew = skew; + P = p; + } + + /// + /// Gets or sets the random number generator which is used to draw random samples. + /// + public System.Random RandomSource + { + get { return _random; } + set { _random = value ?? SystemRandomSource.Default; } + } + + /// + /// A string representation of the distribution. + /// + /// a string representation of the distribution. + public override string ToString() + { + return $"SkewedGeneralizedError(μ = {Location}, σ = {Scale}, λ = { Skew }, p = {P}"; + } + + /// + /// Tests whether the provided values are valid parameters for this distribution. + /// + /// The location (μ) of the distribution. + /// The scale (σ) of the distribution. Range: σ > 0. + /// The skew, 1 > λ > -1 + /// Parameter that controls kurtosis. Range: p > 0 + public static bool IsValidParameterSet(double location, double scale, double skew, double p) + { + return scale > 0.0 && skew > -1.0 && skew < 1.0 && p > 0.0 && !double.IsNaN(location); + } + + /// + /// Gets the location (μ) of the Skewed Generalized t-distribution. + /// + public double Location { get; private set; } + + /// + /// Gets the scale (σ) of the Skewed Generalized t-distribution. Range: σ > 0. + /// + public double Scale { get; private set; } + + /// + /// Gets the skew (λ) of the Skewed Generalized t-distribution. Range: 1 > λ > -1. + /// + public double Skew { get; private set; } + + /// + /// Gets the parameter that controls the kurtosis of the distribution. Range: p > 0. + /// + public double P { get; private set; } + + public double Mode => throw new NotImplementedException(); + + public double Minimum => double.NegativeInfinity; + + public double Maximum => double.PositiveInfinity; + + public double Mean => Location; + + public double Variance => Scale * Scale; + + public double StdDev => Scale; + + public double Entropy => throw new NotImplementedException(); + + public double Skewness => throw new NotImplementedException(); + + public double Median => Location; + + private static double AdjustScale(double scale, double skew, double p) + { + var g1 = SpecialFunctions.Gamma(3.0 / p); + var g2 = SpecialFunctions.Gamma(0.5 + 1.0 / p); + var g3 = SpecialFunctions.Gamma(1.0 / p); + var g4 = SpecialFunctions.Gamma(1.0 / p); + var n1 = Constants.Pi * (1.0 + 3.0 * skew * skew) * g1; + var n2 = Math.Pow(16.0, 1.0 / p) * skew * skew * Math.Pow(g2, 2) * g3; + var d = Constants.Pi * g4; + return scale / Math.Sqrt((n1 - n2) / d); + } + + private static double AdjustX(double x, double scale, double skew, double p) + { + return x + AdjustAddend(scale, skew, p); + } + + private static double AdjustAddend(double scale, double skew, double p) + { + return (Math.Pow(2.0, 2.0 / p) * scale * skew * SpecialFunctions.Gamma(1.0 / 2.0 + 1.0 / p)) / + Math.Sqrt(Constants.Pi); + } + + public static double PDF(double location, double scale, double skew, double p, double x) + { + if (!IsValidParameterSet(location, scale, skew, p)) + { + throw new ArgumentException(Resources.InvalidDistributionParameters); + } + + scale = AdjustScale(scale, skew, p); + x = AdjustX(x, scale, skew, p); + + // p / (2 * sigma * gamma(1 / p) * exp((abs(x - mu) / (sigma * (1 + lambda * sgn(x - mu)))) ^ p)) + var d1 = Math.Abs(x - location); + var d2 = scale * (1.0 + skew * Math.Sign(x - location)); + var d3 = 2.0 * scale * SpecialFunctions.Gamma(1.0 / p); + return p / (Math.Exp(Math.Pow(d1 / d2, p)) * d3); + } + + public static double PDFLn(double location, double scale, double skew, double p, double x) + { + if (!IsValidParameterSet(location, scale, skew, p)) + { + throw new ArgumentException(Resources.InvalidDistributionParameters); + } + + scale = AdjustScale(scale, skew, p); + x = AdjustX(x, scale, skew, p); + + return Math.Log(p) - Math.Log(2.0) - Math.Log(scale) - SpecialFunctions.GammaLn(1.0 / p) - + Math.Pow(Math.Abs(x - location) / (scale * (1.0 + skew * Math.Sign(x - location))), p); + } + + public static double CDF(double location, double scale, double skew, double p, double x) + { + if (!IsValidParameterSet(location, scale, skew, p)) + { + throw new ArgumentException(Resources.InvalidDistributionParameters); + } + + scale = AdjustScale(scale, skew, p); + x = AdjustX(x, scale, skew, p) - location; + + var flip = x < 0; + if (flip) + { + skew = -skew; + x = -x; + } + + var res = (1.0 - skew) / 2.0 + (1.0 + skew) / 2.0 * Gamma.CDF(1.0 / p, 1.0, Math.Pow(x / (scale * (1.0 + skew)), p)); + return flip ? 1.0 - res : res; + } + + public static double InvCDF(double location, double scale, double skew, double p, double pr) + { + if (!IsValidParameterSet(location, scale, skew, p)) + { + throw new ArgumentException(Resources.InvalidDistributionParameters); + } + + scale = AdjustScale(scale, skew, p); + + var flip = pr < (1.0 - skew) / 2.0; + var lambda = skew; + if (flip) + { + pr = 1.0 - pr; + lambda = -lambda; + } + + var res = scale * (1.0 + lambda) * Math.Pow(Gamma.InvCDF(1.0 / p, 1.0, 2 * pr / (1.0 + lambda) + (lambda - 1.0) / (lambda + 1.0)), 1.0 / p); + + if (flip) + res = -res; + res += location; + return res - AdjustAddend(scale, skew, p); + } + + public double InverseCumulativeDistribution(double p) + { + return InvCDF(Location, Scale, Skew, P, p); + } + + public double CumulativeDistribution(double x) + { + return CDF(Location, Scale, Skew, P, x); + } + + public double Density(double x) + { + return PDF(Location, Scale, Skew, P, x); + } + + public double DensityLn(double x) + { + return PDFLn(Location, Scale, Skew, P, x); + } + + /// + /// Generates a sample from the Skew Generalized Error distribution. + /// + /// The random number generator to use. + /// The location (μ) of the distribution. + /// The scale (σ) of the distribution. Range: σ > 0. + /// The skew, 1 > λ > -1 + /// Parameter that controls kurtosis. Range: p > 0 + /// a sample from the distribution. + public static double Sample(System.Random rnd, double location, double scale, double skew, double p) + { + if (!IsValidParameterSet(location, scale, skew, p)) + { + throw new ArgumentException(Resources.InvalidDistributionParameters); + } + + return SampleUnchecked(rnd, location, scale, skew, p); + } + + /// + /// Generates a sequence of samples from the Skew Generalized Error distribution using inverse transform. + /// + /// The random number generator to use. + /// The location (μ) of the distribution. + /// The scale (σ) of the distribution. Range: σ > 0. + /// The skew, 1 > λ > -1 + /// Parameter that controls kurtosis. Range: p > 0 + /// a sequence of samples from the distribution. + public static IEnumerable Samples(System.Random rnd, double location, double scale, double skew, double p) + { + if (!IsValidParameterSet(location, scale, skew, p)) + { + throw new ArgumentException(Resources.InvalidDistributionParameters); + } + + while (true) + { + yield return SampleUnchecked(rnd, location, scale, skew, p); + } + } + + public static IEnumerable Samples(double location, double scale, double skew, double p) + { + if (!IsValidParameterSet(location, scale, skew, p)) + { + throw new ArgumentException(Resources.InvalidDistributionParameters); + } + + return Samples(SystemRandomSource.Default, location, scale, skew, p); + } + + public static double Sample(double location, double scale, double skew, double p) + { + if (!IsValidParameterSet(location, scale, skew, p)) + { + throw new ArgumentException(Resources.InvalidDistributionParameters); + } + + return SampleUnchecked(SystemRandomSource.Default, location, scale, skew, p); + } + + private static double SampleUnchecked(System.Random rnd, double location, double scale, double skew, double p) + { + var u = ContinuousUniform.Sample(rnd, 0, 1); + return InvCDF(location, scale, skew, p, u); + } + + public double Sample() + { + return SampleUnchecked(SystemRandomSource.Default, Location, Scale, Skew, P); + } + + public void Samples(double[] values) + { + if (values == null) + return; + + for (int i = 0; i < values.Length; i++) + { + values[i] = Sample(); + } + } + + public IEnumerable Samples() + { + return Samples(); + } + } +} diff --git a/src/Numerics/Distributions/SkewedGeneralizedT.cs b/src/Numerics/Distributions/SkewedGeneralizedT.cs new file mode 100644 index 00000000..d655221b --- /dev/null +++ b/src/Numerics/Distributions/SkewedGeneralizedT.cs @@ -0,0 +1,523 @@ +// +// 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 System; +using System.Collections.Generic; +using MathNet.Numerics.Properties; +using MathNet.Numerics.Random; + +namespace MathNet.Numerics.Distributions +{ + + /// + /// Continuous Univariate Skewed Generalized T-distribution. + /// Implements the univariate Skewed Generalized t-distribution. For details about this + /// distribution, see + /// + /// Wikipedia - Skewed generalized t-distribution. + /// The skewed generalized t-distribution contains many different distributions within it + /// as special cases based on the parameterization chosen. + /// + /// This implementation is based on the R package dsgt and corresponding viginette, see + /// https://cran.r-project.org/web/packages/sgt/vignettes/sgt.pdf. Compared to that + /// implementation, the options for mean adjustment and variance adjustment are always true. + /// The location (μ) is the mean of the distribution. + /// The scale (σ) squared is the variance of the distribution. + /// + /// The distribution will use the by + /// default. Users can get/set the random number generator by using the + /// property. + /// The statistics classes will check all the incoming parameters + /// whether they are in the allowed range. + public class SkewedGeneralizedT : IContinuousDistribution + { + private System.Random _random; + + // If the given parameterization is one of the recognized special cases, then + // this variable is non-null and the special case is used for all functions. + // Else this value is null and the full formulation of the generalized distribution is used. + private IContinuousDistribution _d; + + /// + /// Initializes a new instance of the SkewedGeneralizedT class. This is a skewed generalized t-distribution + /// with location=0.0, scale=1.0, skew=0.0, p=2.0 and q=Inf (a standard normal distribution). + /// + public SkewedGeneralizedT() + { + _random = SystemRandomSource.Default; + Location = 0.0; + Scale = 1.0; + Skew = 0.0; + P = 2.0; + Q = double.PositiveInfinity; + + _d = new Normal(Location, Scale, _random); + } + + /// + /// Initializes a new instance of the SkewedGeneralizedT class with a particular location, scale, skew + /// and kurtosis parameters. Different parameterizations result in different distributions. + /// + /// The location (μ) of the distribution. + /// The scale (σ) of the distribution. Range: σ > 0. + /// The skew, 1 > λ > -1 + /// First parameter that controls kurtosis. Range: p > 0 + /// Second parameter that controls kurtosis. Range: q > 0 + public SkewedGeneralizedT(double location, double scale, double skew, double p, double q) + { + if (!IsValidParameterSet(location, scale, skew, p, q)) + { + throw new ArgumentException(Resources.InvalidDistributionParameters); + } + + _random = SystemRandomSource.Default; + Location = location; + Scale = scale; + Skew = skew; + P = p; + Q = q; + + _d = FindSpecializedDistribution(location, scale, skew, p, q); + } + + /// + /// Given a parameter set, returns the distribution that matches this parameterization. + /// + /// The location (μ) of the distribution. + /// The scale (σ) of the distribution. Range: σ > 0. + /// The skew, 1 > λ > -1 + /// First parameter that controls kurtosis. Range: p > 0 + /// Second parameter that controls kurtosis. Range: q > 0 + /// Null if no known distribution matches the parameterization, else the distribution. + public static IContinuousDistribution FindSpecializedDistribution(double location, double scale, double skew, double p, double q) + { + if (p == double.PositiveInfinity) + { + scale *= Math.Sqrt(3.0); + return new ContinuousUniform(location - scale, location + scale); + } + + if (q == double.PositiveInfinity) + return new SkewedGeneralizedError(location, scale, skew, p); + + return null; + } + + /// + /// Gets or sets the random number generator which is used to draw random samples. + /// + public System.Random RandomSource + { + get { return _random; } + set { _random = value ?? SystemRandomSource.Default; } + } + + /// + /// A string representation of the distribution. + /// + /// a string representation of the distribution. + public override string ToString() + { + return $"SkewedGeneralizedT(μ = {Location}, σ = {Scale}, λ = { Skew }, p = {P}, q = {Q})"; + } + + /// + /// Tests whether the provided values are valid parameters for this distribution. + /// + /// The location (μ) of the distribution. + /// The scale (σ) of the distribution. Range: σ > 0. + /// The skew, 1 > λ > -1 + /// First parameter that controls kurtosis. Range: p > 0 + /// Second parameter that controls kurtosis. Range: q > 0 + public static bool IsValidParameterSet(double location, double scale, double skew, double p, double q) + { + return scale > 0.0 && skew > -1.0 && skew < 1.0 && p > 0.0 && q > 0.0 && p*q> 2.0 && !double.IsNaN(location); + } + + /// + /// Gets the location (μ) of the Skewed Generalized t-distribution. + /// + public double Location { get; private set; } + + /// + /// Gets the scale (σ) of the Skewed Generalized t-distribution. Range: σ > 0. + /// + public double Scale { get; private set; } + + /// + /// Gets the skew (λ) of the Skewed Generalized t-distribution. Range: 1 > λ > -1. + /// + public double Skew { get; private set; } + + /// + /// Gets the first parameter that controls the kurtosis of the distribution. Range: p > 0. + /// + public double P { get; private set; } + + /// + /// Gets the second parameter that controls the kurtosis of the distribution. Range: q > 0. + /// + public double Q { get; private set; } + + public double Mode => throw new NotImplementedException(); + + public double Minimum => _d == null ? double.NegativeInfinity : _d.Minimum; + + public double Maximum => _d == null ? double.PositiveInfinity : _d.Maximum; + + public double Mean => _d == null ? Location : _d.Mean; + + public double Variance => _d == null ? Scale * Scale : _d.Variance; + + public double StdDev => _d == null ? Scale : _d.StdDev; + + public double Entropy => _d == null ? throw new NotImplementedException() : _d.Entropy; + + public double Skewness => _d == null ? throw new NotImplementedException() : _d.Skewness; + + public double Median => _d == null ? Location : _d.Median; + + private static double AdjustScale(double scale, double skew, double p, double q) + { + var b1 = SpecialFunctions.Beta(3.0 / p, q - 2.0 / p); + var b2 = SpecialFunctions.Beta(1.0 / p, q); + var b3 = SpecialFunctions.Beta(2.0 / p, q - 1.0 / p); + var b4 = SpecialFunctions.Beta(1.0 / p, q); + + return scale / (Math.Pow(q, 1.0 / p) * Math.Sqrt((3.0 * skew * skew + 1.0) * b1 / b2 - 4.0 * skew * skew * ((b3 / b4) * (b3 / b4)))); + } + + // Note: Scale is assumed to be adjusted already when calling this function. + private static double AdjustX(double x, double scale, double skew, double p, double q) + { + return x + AdjustAddend(scale, skew, p, q); + } + + // Note: Scale is assumed to be adjusted already when calling this function. + private static double AdjustAddend(double scale, double skew, double p, double q) + { + var b1 = SpecialFunctions.Beta(2.0 / p, q - 1.0 / p); + var b2 = SpecialFunctions.Beta(1.0 / p, q); + + return (2.0 * scale * skew * Math.Pow(q, 1.0 / p) * b1) / b2; + } + + /// + /// Computes the probability density of the distribution (PDF) at x, i.e. ∂P(X ≤ x)/∂x. + /// + /// The location (μ) of the distribution. + /// The scale (σ) of the distribution. Range: σ > 0. + /// The skew, 1 > λ > -1 + /// First parameter that controls kurtosis. Range: p > 0 + /// Second parameter that controls kurtosis. Range: q > 0 + /// The location at which to compute the density. + /// the density at . + /// + public static double PDF(double location, double scale, double skew, double p, double q, double x) + { + if (!IsValidParameterSet(location, scale, skew, p, q)) + { + throw new ArgumentException(Resources.InvalidDistributionParameters); + } + + var fn = PDFunc(location, scale, skew, p, q, false); + return fn(x); + } + + /// + /// Computes the log probability density of the distribution (lnPDF) at x, i.e. ln(∂P(X ≤ x)/∂x). + /// + /// The location (μ) of the distribution. + /// The scale (σ) of the distribution. Range: σ > 0. + /// The skew, 1 > λ > -1 + /// First parameter that controls kurtosis. Range: p > 0 + /// Second parameter that controls kurtosis. Range: q > 0 + /// The location at which to compute the density. + /// the density at . + /// + public static double PDFLn(double location, double scale, double skew, double p, double q, double x) + { + if (!IsValidParameterSet(location, scale, skew, p, q)) + { + throw new ArgumentException(Resources.InvalidDistributionParameters); + } + + var fn = PDFunc(location, scale, skew, p, q, true); + return fn(x); + } + + private static double PDFull(double location, double scale, double skew, double p, double q, double x) + { + scale = AdjustScale(scale, skew, p, q); + x = AdjustX(x, scale, skew, p, q); + + var b = SpecialFunctions.Beta(1.0 / p, q); + var skewSign = Math.Sign(x - location); + var d1 = Math.Pow(Math.Abs(x - location), p); + var d2 = q * Math.Pow(scale, p) * Math.Pow(skew * skewSign + 1.0, p); + + var denominator = 2.0 * scale * Math.Pow(q, 1.0 / p) * b * Math.Pow(d1 / d2 + 1.0, 1.0 / p + q); + return p / denominator; + } + + private static double PDFullLn(double location, double scale, double skew, double p, double q, double x) + { + scale = AdjustScale(scale, skew, p, q); + x = AdjustX(x, scale, skew, p, q); + + var bLn = SpecialFunctions.BetaLn(1.0 / p, q); + return Math.Log(p) - Math.Log(2.0) - Math.Log(scale) - Math.Log(q) / p - bLn - (1.0 / p + q) * + Math.Log(1.0 + Math.Pow(Math.Abs(x - location), p) / + (q * Math.Pow(scale, p) * Math.Pow(1.0 + skew * Math.Sign(x - location), p))); + } + + // For known parameterizations we just use the existing distributions as visualized + // by Hansen, McDonald and Newey (2010). + // Note that, for all cases where skew is required to be 0, if skew is non-zero, this + // simply gives the corresponding skewed version of the distribution. + private static Func PDFunc(double location, double scale, double skew, double p, double q, bool ln) + { + if (p == double.PositiveInfinity) + { + scale *= Math.Sqrt(3.0); + return x => ln ? ContinuousUniform.PDFLn(location - scale, location + scale, x) : + ContinuousUniform.PDF(-1.0 * (Math.Sqrt(3.0) * scale + location), Math.Sqrt(3.0) * scale + location, x); + } + if (q == double.PositiveInfinity) + return x => ln ? SkewedGeneralizedError.PDFLn(location, scale, skew, p, x) : + SkewedGeneralizedError.PDF(location, scale, skew, p, x); + + return x => ln ? PDFullLn(location, scale, skew, p, q, x) : + PDFull(location, scale, skew, p, q, 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 location (μ) of the distribution. + /// The scale (σ) of the distribution. Range: σ > 0. + /// The skew, 1 > λ > -1 + /// First parameter that controls kurtosis. Range: p > 0 + /// Second parameter that controls kurtosis. Range: q > 0 + /// the cumulative distribution at location . + /// + public static double CDF(double location, double scale, double skew, double p, double q, double x) + { + if (!IsValidParameterSet(location, scale, skew, p, q)) + { + throw new ArgumentException(Resources.InvalidDistributionParameters); + } + + // Note: Adapted from the R package, + // based on a transformation of the cumulative probability density function that uses the + // incomplete beta function or incomplete gamma function. + + scale = AdjustScale(scale, skew, p, q); + x = AdjustX(x, scale, skew, p, q) - location; + + var flip = x > 0; + if (flip) + { + skew = -skew; + x = -x; + } + + var res = (1.0 - skew) / 2.0 + (skew - 1.0) / 2.0 * Beta.CDF(1.0 / p, q, 1.0 / (1.0 + q * Math.Pow(scale * (1.0 - skew) / -x, p))); + return flip ? 1.0 - res : res; + } + + /// + /// 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 (μ) of the distribution. + /// The scale (σ) of the distribution. Range: σ > 0. + /// The skew, 1 > λ > -1 + /// First parameter that controls kurtosis. Range: p > 0 + /// Second parameter that controls kurtosis. Range: q > 0 + /// the inverse cumulative density at . + /// + public static double InvCDF(double location, double scale, double skew, double p, double q, double pr) + { + if (!IsValidParameterSet(location, scale, skew, p, q)) + { + throw new ArgumentException(Resources.InvalidDistributionParameters); + } + + // Note: Adapted from the R package, + // solving for the inverse of the CDF that uses the inverse of the incomplete beta function or + // incomplete gamma function + + scale = AdjustScale(scale, skew, p, q); + + var flip = pr > (1.0 - skew) / 2.0; + var lambda = skew; + if (flip) + { + pr = 1.0 - pr; + lambda = -lambda; + } + + var res = scale * (lambda - 1.0) * Math.Pow(1.0 / (q * Beta.InvCDF(1.0 / p, q, 1.0 - 2.0 * pr / (1.0 - lambda))) - 1.0 / q, -1.0 / p); + + if (flip) + res = -res; + res += location; + return res - AdjustAddend(scale, skew, p, q); + } + + public double CumulativeDistribution(double x) + { + return _d == null ? CDF(Location, Scale, Skew, P, Q, x) : _d.CumulativeDistribution(x); + } + + /// + /// Computes the inverse of the cumulative distribution function (InvCDF) for the distribution + /// at the given probability. This is also known as the quantile or percent point function. + /// + /// The location at which to compute the inverse cumulative density. + /// the inverse cumulative density at . + /// + public double InverseCumulativeDistribution(double p) + { + // InverseCumulativeDistribution is not a part of the interface, so resort to type-checking. + if (_d != null) + { + if (_d is SkewedGeneralizedError sge) + return sge.InverseCumulativeDistribution(p); + if (_d is ContinuousUniform u) + return u.InverseCumulativeDistribution(p); + } + + return InvCDF(Location, Scale, Skew, P, Q, p); + } + + public double Density(double x) + { + return _d == null ? PDF(Location, Scale, Skew, P, Q, x) : _d.Density(x); + } + + public double DensityLn(double x) + { + return _d == null ? PDFLn(Location, Scale, Skew, P, Q, x) : _d.DensityLn(x); + } + + /// + /// Generates a sample from the Skew Generalized t-distribution. + /// + /// The random number generator to use. + /// The location (μ) of the distribution. + /// The scale (σ) of the distribution. Range: σ > 0. + /// The skew, 1 > λ > -1 + /// First parameter that controls kurtosis. Range: p > 0 + /// Second parameter that controls kurtosis. Range: q > 0 + /// a sample from the distribution. + public static double Sample(System.Random rnd, double location, double scale, double skew, double p, double q) + { + if (!IsValidParameterSet(location, scale, skew, p, q)) + { + throw new ArgumentException(Resources.InvalidDistributionParameters); + } + + return SampleUnchecked(rnd, location, scale, skew, p, q); + } + + /// + /// Generates a sequence of samples from the Skew Generalized t-distribution using inverse transform. + /// + /// The random number generator to use. + /// The location (μ) of the distribution. + /// The scale (σ) of the distribution. Range: σ > 0. + /// The skew, 1 > λ > -1 + /// First parameter that controls kurtosis. Range: p > 0 + /// Second parameter that controls kurtosis. Range: q > 0 + /// a sequence of samples from the distribution. + public static IEnumerable Samples(System.Random rnd, double location, double scale, double skew, double p, double q) + { + if (!IsValidParameterSet(location, scale, skew, p, q)) + { + throw new ArgumentException(Resources.InvalidDistributionParameters); + } + + while (true) + { + yield return SampleUnchecked(rnd, location, scale, skew, p, q); + } + } + + public static IEnumerable Samples(double location, double scale, double skew, double p, double q) + { + if (!IsValidParameterSet(location, scale, skew, p, q)) + { + throw new ArgumentException(Resources.InvalidDistributionParameters); + } + + return Samples(SystemRandomSource.Default, location, scale, skew, p, q); + } + + public static double Sample(double location, double scale, double skew, double p, double q) + { + if (!IsValidParameterSet(location, scale, skew, p, q)) + { + throw new ArgumentException(Resources.InvalidDistributionParameters); + } + + return SampleUnchecked(SystemRandomSource.Default, location, scale, skew, p, q); + } + + private static double SampleUnchecked(System.Random rnd, double location, double scale, double skew, double p, double q) + { + var u = ContinuousUniform.Sample(rnd, 0, 1); + return InvCDF(location, scale, skew, p, q, u); + } + + public double Sample() + { + return SampleUnchecked(SystemRandomSource.Default, Location, Scale, Skew, P, Q); + } + + public void Samples(double[] values) + { + if (values == null) + return; + + for (int i = 0; i < values.Length; i++) + { + values[i] = Sample(); + } + } + + public IEnumerable Samples() + { + return Samples(); + } + } +}