From 04a8fa8dc5e7fff2a40945f70afd59f66713585e Mon Sep 17 00:00:00 2001 From: mikael Date: Sun, 11 Aug 2019 18:19:26 +0200 Subject: [PATCH] Add skewness and mode calculation --- .../Continuous/SkewGeneralizedTTests.cs | 62 +++++++++++++++++++ .../Distributions/SkewedGeneralizedError.cs | 36 ++++++++++- .../Distributions/SkewedGeneralizedT.cs | 42 ++++++++++++- 3 files changed, 134 insertions(+), 6 deletions(-) diff --git a/src/Numerics.Tests/DistributionTests/Continuous/SkewGeneralizedTTests.cs b/src/Numerics.Tests/DistributionTests/Continuous/SkewGeneralizedTTests.cs index 659e43da..f656afba 100644 --- a/src/Numerics.Tests/DistributionTests/Continuous/SkewGeneralizedTTests.cs +++ b/src/Numerics.Tests/DistributionTests/Continuous/SkewGeneralizedTTests.cs @@ -218,6 +218,68 @@ namespace MathNet.Numerics.UnitTests.DistributionTests.Continuous Assert.IsTrue(sp < p); } + [TestCase(0, 1, -0.1, 0.5123)] + [TestCase(0, 1, 0.1, 0.6123)] + public void ValidateModeOfSkewedNormalDistribution(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 sm = sn.Mode; + var m = n.Mode; + + if (skew < 0) + Assert.IsTrue(sm > m); + else + Assert.IsTrue(sm < m); + } + + [TestCase(0, 1, -0.1, 3, 5, 0.5123)] + [TestCase(0, 1, 0.1, 3, 5, 0.6123)] + public void ValidateModeOfSkewedGeneralizedTDistribution(double location, double scale, double skew, double p, double q, double x) + { + var sn = new SkewedGeneralizedT(location, scale, skew, 2, double.PositiveInfinity); + var n = new Normal(location, scale); + + var sm = sn.Mode; + var m = n.Mode; + + if (skew < 0) + Assert.IsTrue(sm > m); + else + Assert.IsTrue(sm < m); + } + + [TestCase(-0.5)] + [TestCase(0.0)] + [TestCase(0.5)] + public void ValidateSkewnessOfNormalDistribution(double skew) + { + var sn = new SkewedGeneralizedT(0.0, 1.0, skew, 2.0, double.PositiveInfinity); + + if (skew > 0) + Assert.IsTrue(sn.Skewness > 0.0); + else if (skew < 0) + Assert.IsTrue(sn.Skewness < 0.0); + else + Assert.AreEqual(skew, sn.Skewness); + } + + [TestCase(-0.5)] + [TestCase(0.0)] + [TestCase(0.5)] + public void ValidateSkewnessOfSkewGeneralizedTDistribution(double skew) + { + var sn = new SkewedGeneralizedT(0.0, 1.0, skew, 3.0, 4.0); + + if (skew > 0) + Assert.IsTrue(sn.Skewness > 0.0); + else if (skew < 0) + Assert.IsTrue(sn.Skewness < 0.0); + else + Assert.AreEqual(skew, sn.Skewness); + } + /// /// Can sample static. /// diff --git a/src/Numerics/Distributions/SkewedGeneralizedError.cs b/src/Numerics/Distributions/SkewedGeneralizedError.cs index 78f57fd7..baeb5167 100644 --- a/src/Numerics/Distributions/SkewedGeneralizedError.cs +++ b/src/Numerics/Distributions/SkewedGeneralizedError.cs @@ -58,6 +58,8 @@ namespace MathNet.Numerics.Distributions { private System.Random _random; + private readonly double _skewness; + /// /// 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). @@ -91,6 +93,8 @@ namespace MathNet.Numerics.Distributions Scale = scale; Skew = skew; P = p; + + _skewness = CalculateSkewness(); } /// @@ -143,23 +147,49 @@ namespace MathNet.Numerics.Distributions /// public double P { get; private set; } - public double Mode => throw new NotImplementedException(); + // No skew implies Median=Mode=Mean + public double Mode => + Skew == 0 ? Mean : Mean - AdjustAddend(AdjustScale(Scale, Skew, P), Skew, P); public double Minimum => double.NegativeInfinity; public double Maximum => double.PositiveInfinity; + // Mean=Location due to our adjustments made public double Mean => Location; + // Variance=Scale*Scale due to our adjustments made public double Variance => Scale * Scale; public double StdDev => Scale; public double Entropy => throw new NotImplementedException(); - public double Skewness => throw new NotImplementedException(); + public double Skewness => _skewness; + + // No skew implies Median=Mode=Mean + // Else find it via the point where CDF gives 0.5 + public double Median => + Skew == 0 ? Mean : InverseCumulativeDistribution(0.5); + + private double CalculateSkewness() + { + if (Skew == 0) + return 0.0; + + var piPow = Math.Pow(Constants.Pi, 3.0 / 2.0); + var g1 = SpecialFunctions.Gamma(1.0 / P); + var g2 = SpecialFunctions.Gamma(0.5 + 1.0 / P); + var g3 = SpecialFunctions.Gamma(3.0 / P); + var g4 = SpecialFunctions.Gamma(4.0 / P); - public double Median => Location; + var t1 = Skew * Scale * Scale * Scale / (piPow * g1); + var t2 = Math.Pow(2.0, (6.0 + P) / P) * Skew * Skew * Math.Pow(g2, 3.0) * g1; + var t3 = 3.0 * Math.Pow(4.0, 1.0 / P) * Constants.Pi * (1.0 + 3.0 * Skew * Skew) * g2 * g3; + var t4 = 4.0 * piPow * (1.0 + Skew * Skew) * g4; + + return t1 * (t2 - t3 + t4); + } private static double AdjustScale(double scale, double skew, double p) { diff --git a/src/Numerics/Distributions/SkewedGeneralizedT.cs b/src/Numerics/Distributions/SkewedGeneralizedT.cs index d655221b..720cdeec 100644 --- a/src/Numerics/Distributions/SkewedGeneralizedT.cs +++ b/src/Numerics/Distributions/SkewedGeneralizedT.cs @@ -64,6 +64,8 @@ namespace MathNet.Numerics.Distributions // Else this value is null and the full formulation of the generalized distribution is used. private IContinuousDistribution _d; + private readonly double _skewness; + /// /// 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). @@ -104,6 +106,9 @@ namespace MathNet.Numerics.Distributions Q = q; _d = FindSpecializedDistribution(location, scale, skew, p, q); + + if (_d == null) + _skewness = CalculateSkewness(); } /// @@ -185,23 +190,54 @@ namespace MathNet.Numerics.Distributions /// public double Q { get; private set; } - public double Mode => throw new NotImplementedException(); + // No skew implies Median=Mode=Mean + public double Mode => _d == null ? + Skew == 0 ? Mean : Mean - AdjustAddend(AdjustScale(Scale, Skew, P, Q), Skew, P, Q) : + _d.Mode; public double Minimum => _d == null ? double.NegativeInfinity : _d.Minimum; public double Maximum => _d == null ? double.PositiveInfinity : _d.Maximum; + // Mean=Location due to our adjustments made public double Mean => _d == null ? Location : _d.Mean; + // Variance=Scale*Scale due to our adjustments made 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 Skewness => _d == null ? + _skewness : + _d.Skewness; + + // No skew implies Median=Mode=Mean + // Else find it via the point where CDF gives 0.5 + public double Median => _d == null ? + Skew == 0 ? Mean : InverseCumulativeDistribution(0.5) : + _d.Median; - public double Median => _d == null ? Location : _d.Median; + private double CalculateSkewness() + { + if (P * Q <= 3 || Skew == 0) + return 0.0; + + var scale = AdjustScale(Scale, Skew, P, Q); + var b1 = SpecialFunctions.Beta(1.0 / P, Q); + var b2 = SpecialFunctions.Beta(2.0 / P, Q - 1.0 / P); + var b3 = SpecialFunctions.Beta(3.0 / P, Q - 2.0 / P); + var b4 = SpecialFunctions.Beta(4.0 / P, Q - 3.0 / P); + + var t1 = (2.0 * Math.Pow(Q, 3.0 / P) * Skew * Math.Pow(scale, 3.0)) / Math.Pow(b1, 3.0); + var t2 = 8.0 * Skew * Skew * Math.Pow(b2, 3.0); + var t3 = 3.0 * (1.0 + 3.0 * Skew * Skew) * b1; + var t4 = b2 * b3; + var t5 = 2.0 * (1.0 + Skew * Skew) * Math.Pow(b1, 2.0) * b4; + + return t1 * (t2 - t3 * t4 + t5); + } private static double AdjustScale(double scale, double skew, double p, double q) {