Add skewness and mode calculation

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);
+        }
+
         /// <summary>
         /// Can sample static.
         /// </summary>

src/Numerics/Distributions/SkewedGeneralizedError.cs

@@ -58,6 +58,8 @@ namespace MathNet.Numerics.Distributions
     {
         private System.Random _random;
 
+        private readonly double _skewness;
+
         /// <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).
@@ -91,6 +93,8 @@ namespace MathNet.Numerics.Distributions
             Scale = scale;
             Skew = skew;
             P = p;
+
+            _skewness = CalculateSkewness();
         }
 
         /// <summary>
@@ -143,23 +147,49 @@ namespace MathNet.Numerics.Distributions
         /// </summary>
         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)
         {

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;
+
         /// <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).
@@ -104,6 +106,9 @@ namespace MathNet.Numerics.Distributions
             Q = q;
 
             _d = FindSpecializedDistribution(location, scale, skew, p, q);
+
+            if (_d == null)
+                _skewness = CalculateSkewness();
         }
 
         /// <summary>
@@ -185,23 +190,54 @@ namespace MathNet.Numerics.Distributions
         /// </summary>
         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)
         {