uo9xpb74
/
mathnet-numerics
forked from tsai/mathnet-numerics
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity
Browse Source

Add Skewed Generalized T distribution and Skewed Generalized Error distribution

v4
mikael 7 years ago
committed by Christoph Ruegg
parent
982966f704
commit
332251d99c
4 changed files with 1164 additions and 0 deletions
Whitespace
Show all changes Ignore whitespace when comparing lines Ignore changes in amount of whitespace Ignore changes in whitespace at EOL
Unified View
Diff Options
Show Stats Download Patch File Download Diff File
  1. 8
      src/FSharp/Distributions.fs
  2. 261
      src/Numerics.Tests/DistributionTests/Continuous/SkewGeneralizedTTests.cs
  3. 372
      src/Numerics/Distributions/SkewedGeneralizedError.cs
  4. 523
      src/Numerics/Distributions/SkewedGeneralizedT.cs

8
src/FSharp/Distributions.fs
View File

@ -150,6 +150,14 @@ module Sample =
let studentT location scale freedom (rng:System.Random) = StudentT.Sample(rng, location, scale, freedom) 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) 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 (λ). /// Weibull with shape (k) and scale (λ).
let weibull shape scale (rng:System.Random) = Weibull.Sample(rng, shape, 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) let weibullSeq shape scale (rng:System.Random) = Weibull.Samples(rng, shape, scale)

261
src/Numerics.Tests/DistributionTests/Continuous/SkewGeneralizedTTests.cs
View File

@ -0,0 +1,261 @@
// <copyright file="StudentTTests.cs" company="Math.NET">
// 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.
// </copyright>
using System;
using System.Linq;
using MathNet.Numerics.Distributions;
using NUnit.Framework;
namespace MathNet.Numerics.UnitTests.DistributionTests.Continuous
{
/// <summary>
/// <c>SkewedGeneralizedT</c> distribution tests.
/// Reference values are from the R package sgt 2.0 (run on Microsoft R Open v3.5.2)
/// </summary>
[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);
}
/// <summary>
/// Can sample static.
/// </summary>
[Test]
public void CanSampleStatic()
{
SkewedGeneralizedT.Sample(0.0, 1.0, 0.3, 2.2, 5.6);
}
/// <summary>
/// Can sample sequence static.
/// </summary>
[Test]
public void CanSampleSequenceStatic()
{
var ied = SkewedGeneralizedT.Samples(0.0, 1.0, 0.3, 2.2, 5.6);
GC.KeepAlive(ied.Take(5).ToArray());
}
/// <summary>
/// Can sample.
/// </summary>
[Test]
public void CanSample()
{
var n = new SkewedGeneralizedT();
n.Sample();
}
/// <summary>
/// Can sample sequence.
/// </summary>
[Test]
public void CanSampleSequence()
{
var n = new SkewedGeneralizedT();
var ied = n.Samples();
GC.KeepAlive(ied.Take(5).ToArray());
}
}
}

372
src/Numerics/Distributions/SkewedGeneralizedError.cs
View File

@ -0,0 +1,372 @@
// <copyright file="StudentT.cs" company="Math.NET">
// 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.
// </copyright>
using MathNet.Numerics.Properties;
using MathNet.Numerics.Random;
using System;
using System.Collections.Generic;
namespace MathNet.Numerics.Distributions
{
/// <summary>
/// Continuous Univariate Skewed Generalized Error Distribution (SGED).
/// Implements the univariate SSkewed Generalized Error Distribution. For details about this
/// distribution, see
/// <a href="https://en.wikipedia.org/wiki/Generalized_normal_distribution">
/// Wikipedia - Generalized Error Distribution</a>.
/// It includes Laplace, Normal and Student-t distributions.
/// This is the <see cref="SkewedGeneralizedT"/> distribution with q=Inf.
/// </summary>
/// <remarks><para>This implementation is based on the R package dsgt and corresponding viginette, see
/// <a href="">https://cran.r-project.org/web/packages/sgt/vignettes/sgt.pdf</a>. 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.
/// </para>
/// <para>The distribution will use the <see cref="System.Random"/> by
/// default. Users can get/set the random number generator by using the
/// <see cref="RandomSource"/> property.</para>
/// <para>The statistics classes will check all the incoming parameters
/// whether they are in the allowed range.</para></remarks>
public class SkewedGeneralizedError : IContinuousDistribution
{
private System.Random _random;
/// <summary>
/// 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).
/// </summary>
public SkewedGeneralizedError()
{
_random = SystemRandomSource.Default;
Location = 0.0;
Scale = 1.0;
Skew = 0.0;
P = 2.0;
}
/// <summary>
/// Initializes a new instance of the SkewedGeneralizedT class with a particular location, scale, skew
/// and kurtosis parameters. Different parameterizations result in different distributions.
/// </summary>
/// <param name="location">The location (μ) of the distribution.</param>
/// <param name="scale">The scale (σ) of the distribution. Range: σ > 0.</param>
/// <param name="skew">The skew, 1 > λ > -1</param>
/// <param name="p">Parameter that controls kurtosis. Range: p > 0</param>
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;
}
/// <summary>
/// Gets or sets the random number generator which is used to draw random samples.
/// </summary>
public System.Random RandomSource
{
get { return _random; }
set { _random = value ?? SystemRandomSource.Default; }
}
/// <summary>
/// A string representation of the distribution.
/// </summary>
/// <returns>a string representation of the distribution.</returns>
public override string ToString()
{
return $"SkewedGeneralizedError(μ = {Location}, σ = {Scale}, λ = { Skew }, p = {P}";
}
/// <summary>
/// Tests whether the provided values are valid parameters for this distribution.
/// </summary>
/// <param name="location">The location (μ) of the distribution.</param>
/// <param name="scale">The scale (σ) of the distribution. Range: σ > 0.</param>
/// <param name="skew">The skew, 1 > λ > -1</param>
/// <param name="p">Parameter that controls kurtosis. Range: p > 0</param>
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);
}
/// <summary>
/// Gets the location (μ) of the Skewed Generalized t-distribution.
/// </summary>
public double Location { get; private set; }
/// <summary>
/// Gets the scale (σ) of the Skewed Generalized t-distribution. Range: σ > 0.
/// </summary>
public double Scale { get; private set; }
/// <summary>
/// Gets the skew (λ) of the Skewed Generalized t-distribution. Range: 1 > λ > -1.
/// </summary>
public double Skew { get; private set; }
/// <summary>
/// Gets the parameter that controls the kurtosis of the distribution. Range: p > 0.
/// </summary>
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);
}
/// <summary>
/// Generates a sample from the Skew Generalized Error distribution.
/// </summary>
/// <param name="rnd">The random number generator to use.</param>
/// <param name="location">The location (μ) of the distribution.</param>
/// <param name="scale">The scale (σ) of the distribution. Range: σ > 0.</param>
/// <param name="skew">The skew, 1 > λ > -1</param>
/// <param name="p">Parameter that controls kurtosis. Range: p > 0</param>
/// <returns>a sample from the distribution.</returns>
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);
}
/// <summary>
/// Generates a sequence of samples from the Skew Generalized Error distribution using inverse transform.
/// </summary>
/// <param name="rnd">The random number generator to use.</param>
/// <param name="location">The location (μ) of the distribution.</param>
/// <param name="scale">The scale (σ) of the distribution. Range: σ > 0.</param>
/// <param name="skew">The skew, 1 > λ > -1</param>
/// <param name="p">Parameter that controls kurtosis. Range: p > 0</param>
/// <returns>a sequence of samples from the distribution.</returns>
public static IEnumerable<double> 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<double> 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<double> Samples()
{
return Samples();
}
}
}

523
src/Numerics/Distributions/SkewedGeneralizedT.cs
View File

@ -0,0 +1,523 @@
// <copyright file="StudentT.cs" company="Math.NET">
// 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.
// </copyright>
using System;
using System.Collections.Generic;
using MathNet.Numerics.Properties;
using MathNet.Numerics.Random;
namespace MathNet.Numerics.Distributions
{
/// <summary>
/// Continuous Univariate Skewed Generalized T-distribution.
/// Implements the univariate Skewed Generalized t-distribution. For details about this
/// distribution, see
/// <a href="https://en.wikipedia.org/wiki/Skewed_generalized_t_distribution">
/// Wikipedia - Skewed generalized t-distribution</a>.
/// The skewed generalized t-distribution contains many different distributions within it
/// as special cases based on the parameterization chosen.
/// </summary>
/// <remarks><para>This implementation is based on the R package dsgt and corresponding viginette, see
/// <a href="">https://cran.r-project.org/web/packages/sgt/vignettes/sgt.pdf</a>. 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.
/// </para>
/// <para>The distribution will use the <see cref="System.Random"/> by
/// default. Users can get/set the random number generator by using the
/// <see cref="RandomSource"/> property.</para>
/// <para>The statistics classes will check all the incoming parameters
/// whether they are in the allowed range.</para></remarks>
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;
/// <summary>
/// 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).
/// </summary>
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);
}
/// <summary>
/// Initializes a new instance of the SkewedGeneralizedT class with a particular location, scale, skew
/// and kurtosis parameters. Different parameterizations result in different distributions.
/// </summary>
/// <param name="location">The location (μ) of the distribution.</param>
/// <param name="scale">The scale (σ) of the distribution. Range: σ > 0.</param>
/// <param name="skew">The skew, 1 > λ > -1</param>
/// <param name="p">First parameter that controls kurtosis. Range: p > 0</param>
/// <param name="q">Second parameter that controls kurtosis. Range: q > 0</param>
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);
}
/// <summary>
/// Given a parameter set, returns the distribution that matches this parameterization.
/// </summary>
/// <param name="location">The location (μ) of the distribution.</param>
/// <param name="scale">The scale (σ) of the distribution. Range: σ > 0.</param>
/// <param name="skew">The skew, 1 > λ > -1</param>
/// <param name="p">First parameter that controls kurtosis. Range: p > 0</param>
/// <param name="q">Second parameter that controls kurtosis. Range: q > 0</param>
/// <returns>Null if no known distribution matches the parameterization, else the distribution.</returns>
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;
}
/// <summary>
/// Gets or sets the random number generator which is used to draw random samples.
/// </summary>
public System.Random RandomSource
{
get { return _random; }
set { _random = value ?? SystemRandomSource.Default; }
}
/// <summary>
/// A string representation of the distribution.
/// </summary>
/// <returns>a string representation of the distribution.</returns>
public override string ToString()
{
return $"SkewedGeneralizedT(μ = {Location}, σ = {Scale}, λ = { Skew }, p = {P}, q = {Q})";
}
/// <summary>
/// Tests whether the provided values are valid parameters for this distribution.
/// </summary>
/// <param name="location">The location (μ) of the distribution.</param>
/// <param name="scale">The scale (σ) of the distribution. Range: σ > 0.</param>
/// <param name="skew">The skew, 1 > λ > -1</param>
/// <param name="p">First parameter that controls kurtosis. Range: p > 0</param>
/// <param name="q">Second parameter that controls kurtosis. Range: q > 0</param>
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);
}
/// <summary>
/// Gets the location (μ) of the Skewed Generalized t-distribution.
/// </summary>
public double Location { get; private set; }
/// <summary>
/// Gets the scale (σ) of the Skewed Generalized t-distribution. Range: σ > 0.
/// </summary>
public double Scale { get; private set; }
/// <summary>
/// Gets the skew (λ) of the Skewed Generalized t-distribution. Range: 1 > λ > -1.
/// </summary>
public double Skew { get; private set; }
/// <summary>
/// Gets the first parameter that controls the kurtosis of the distribution. Range: p > 0.
/// </summary>
public double P { get; private set; }
/// <summary>
/// Gets the second parameter that controls the kurtosis of the distribution. Range: q > 0.
/// </summary>
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;
}
/// <summary>
/// Computes the probability density of the distribution (PDF) at x, i.e. ∂P(X ≤ x)/∂x.
/// </summary>
/// <param name="location">The location (μ) of the distribution.</param>
/// <param name="scale">The scale (σ) of the distribution. Range: σ > 0.</param>
/// <param name="skew">The skew, 1 > λ > -1</param>
/// <param name="p">First parameter that controls kurtosis. Range: p > 0</param>
/// <param name="q">Second parameter that controls kurtosis. Range: q > 0</param>
/// <param name="x">The location at which to compute the density.</param>
/// <returns>the density at <paramref name="x"/>.</returns>
/// <seealso cref="Density"/>
public static double PDF(double 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);
}
/// <summary>
/// Computes the log probability density of the distribution (lnPDF) at x, i.e. ln(∂P(X ≤ x)/∂x).
/// </summary>
/// <param name="location">The location (μ) of the distribution.</param>
/// <param name="scale">The scale (σ) of the distribution. Range: σ > 0.</param>
/// <param name="skew">The skew, 1 > λ > -1</param>
/// <param name="p">First parameter that controls kurtosis. Range: p > 0</param>
/// <param name="q">Second parameter that controls kurtosis. Range: q > 0</param>
/// <param name="x">The location at which to compute the density.</param>
/// <returns>the density at <paramref name="x"/>.</returns>
/// <seealso cref="Density"/>
public static double 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<double, double> 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);
}
/// <summary>
/// Computes the cumulative distribution (CDF) of the distribution at x, i.e. P(X ≤ x).
/// </summary>
/// <param name="x">The location at which to compute the cumulative distribution function.</param>
/// <param name="location">The location (μ) of the distribution.</param>
/// <param name="scale">The scale (σ) of the distribution. Range: σ > 0.</param>
/// <param name="skew">The skew, 1 > λ > -1</param>
/// <param name="p">First parameter that controls kurtosis. Range: p > 0</param>
/// <param name="q">Second parameter that controls kurtosis. Range: q > 0</param>
/// <returns>the cumulative distribution at location <paramref name="x"/>.</returns>
/// <seealso cref="CumulativeDistribution"/>
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;
}
/// <summary>
/// Computes the inverse of the cumulative distribution function (InvCDF) for the distribution
/// at the given probability. This is also known as the quantile or percent point function.
/// </summary>
/// <param name="pr">The location at which to compute the inverse cumulative density.</param>
/// <param name="location">The location (μ) of the distribution.</param>
/// <param name="scale">The scale (σ) of the distribution. Range: σ > 0.</param>
/// <param name="skew">The skew, 1 > λ > -1</param>
/// <param name="p">First parameter that controls kurtosis. Range: p > 0</param>
/// <param name="q">Second parameter that controls kurtosis. Range: q > 0</param>
/// <returns>the inverse cumulative density at <paramref name="p"/>.</returns>
/// <seealso cref="InverseCumulativeDistribution"/>
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);
}
/// <summary>
/// Computes the inverse of the cumulative distribution function (InvCDF) for the distribution
/// at the given probability. This is also known as the quantile or percent point function.
/// </summary>
/// <param name="p">The location at which to compute the inverse cumulative density.</param>
/// <returns>the inverse cumulative density at <paramref name="p"/>.</returns>
/// <seealso cref="InvCDF"/>
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);
}
/// <summary>
/// Generates a sample from the Skew Generalized t-distribution.
/// </summary>
/// <param name="rnd">The random number generator to use.</param>
/// <param name="location">The location (μ) of the distribution.</param>
/// <param name="scale">The scale (σ) of the distribution. Range: σ > 0.</param>
/// <param name="skew">The skew, 1 > λ > -1</param>
/// <param name="p">First parameter that controls kurtosis. Range: p > 0</param>
/// <param name="q">Second parameter that controls kurtosis. Range: q > 0</param>
/// <returns>a sample from the distribution.</returns>
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);
}
/// <summary>
/// Generates a sequence of samples from the Skew Generalized t-distribution using inverse transform.
/// </summary>
/// <param name="rnd">The random number generator to use.</param>
/// <param name="location">The location (μ) of the distribution.</param>
/// <param name="scale">The scale (σ) of the distribution. Range: σ > 0.</param>
/// <param name="skew">The skew, 1 > λ > -1</param>
/// <param name="p">First parameter that controls kurtosis. Range: p > 0</param>
/// <param name="q">Second parameter that controls kurtosis. Range: q > 0</param>
/// <returns>a sequence of samples from the distribution.</returns>
public static IEnumerable<double> 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<double> 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<double> Samples()
{
return Samples();
}
}
}
Write Preview
Loading…
Cancel
Save