From 8f91ec929f16af3e53204efbceba6e46df86b47b Mon Sep 17 00:00:00 2001 From: Christoph Ruegg Date: Sat, 4 Jan 2014 20:03:42 +0100 Subject: [PATCH] Distributions: StudentT.InvCDF --- src/Numerics/Distributions/StudentT.cs | 43 +++++++++++++++++++ .../Continuous/NormalTests.cs | 16 +++---- .../Continuous/StudentTTests.cs | 35 ++++++++++----- 3 files changed, 75 insertions(+), 19 deletions(-) diff --git a/src/Numerics/Distributions/StudentT.cs b/src/Numerics/Distributions/StudentT.cs index c52f4e36..4d4b50ce 100644 --- a/src/Numerics/Distributions/StudentT.cs +++ b/src/Numerics/Distributions/StudentT.cs @@ -32,6 +32,7 @@ using System; using System.Collections.Generic; using MathNet.Numerics.Properties; using MathNet.Numerics.Random; +using MathNet.Numerics.RootFinding; namespace MathNet.Numerics.Distributions { @@ -304,6 +305,19 @@ namespace MathNet.Numerics.Distributions return CDF(_location, _scale, _freedom, 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 . + /// + /// WARNING: currently not an explicit implementation, hence slow and unreliable. + public double InverseCumulativeDistribution(double p) + { + return InvCDF(_location, _scale, _freedom, p); + } + /// /// Samples student-t distributed random variables. /// @@ -409,6 +423,35 @@ namespace MathNet.Numerics.Distributions return x <= location ? ib : 1.0 - ib; } + /// + /// 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 degrees of freedom (ν) for the distribution. Range: ν > 0. + /// the inverse cumulative density at . + /// + /// WARNING: currently not an explicit implementation, hence slow and unreliable. + public static double InvCDF(double location, double scale, double freedom, double p) + { + if (scale <= 0.0 || freedom <= 0.0) throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters); + + // TODO JVG we can probably do a better job for Cauchy special case + if (Double.IsPositiveInfinity(freedom)) return Normal.InvCDF(location, scale, p); + if (p == 0.5d) return location; + + // TODO PERF: We must implement this explicitly instead of solving for CDF^-1 + return Brent.FindRoot(x => + { + var k = (x - location)/scale; + var h = freedom/(freedom + (k*k)); + var ib = 0.5*SpecialFunctions.BetaRegularized(freedom/2.0, 0.5, h); + return x <= location ? ib - p : 1.0 - ib - p; + }, -800, 800, accuracy: 1e-8); + } + /// /// Generates a sample from the Student t-distribution. /// diff --git a/src/UnitTests/DistributionTests/Continuous/NormalTests.cs b/src/UnitTests/DistributionTests/Continuous/NormalTests.cs index 30e4d57f..bf163821 100644 --- a/src/UnitTests/DistributionTests/Continuous/NormalTests.cs +++ b/src/UnitTests/DistributionTests/Continuous/NormalTests.cs @@ -458,7 +458,7 @@ namespace MathNet.Numerics.UnitTests.DistributionTests.Continuous /// Validate cumulative distribution. /// /// Input X value. - /// Expected value. + /// Expected value. [TestCase(Double.NegativeInfinity, 0.0)] [TestCase(-5.0, 0.00000028665157187919391167375233287464535385442301361187883)] [TestCase(-2.0, 0.0002326290790355250363499258867279847735487493358890356)] @@ -469,18 +469,18 @@ namespace MathNet.Numerics.UnitTests.DistributionTests.Continuous [TestCase(6.0, 0.69146246127401310363770461060833773988360217555457859)] [TestCase(10.0, 0.9937903346742238648330218954258077788721022530769078)] [TestCase(Double.PositiveInfinity, 1.0)] - public void ValidateCumulativeDistribution(double x, double f) + public void ValidateCumulativeDistribution(double x, double p) { var n = Normal.WithMeanStdDev(5.0, 2.0); - AssertHelpers.AlmostEqualRelative(f, n.CumulativeDistribution(x), 9); - AssertHelpers.AlmostEqualRelative(f, Normal.CDF(5.0, 2.0, x), 9); + AssertHelpers.AlmostEqualRelative(p, n.CumulativeDistribution(x), 9); + AssertHelpers.AlmostEqualRelative(p, Normal.CDF(5.0, 2.0, x), 9); } /// /// Validate inverse cumulative distribution. /// /// Input X value. - /// Expected value. + /// Expected value. [TestCase(Double.NegativeInfinity, 0.0)] [TestCase(-5.0, 0.00000028665157187919391167375233287464535385442301361187883)] [TestCase(-2.0, 0.0002326290790355250363499258867279847735487493358890356)] @@ -491,11 +491,11 @@ namespace MathNet.Numerics.UnitTests.DistributionTests.Continuous [TestCase(6.0, .69146246127401310363770461060833773988360217555457859)] [TestCase(10.0, 0.9937903346742238648330218954258077788721022530769078)] [TestCase(Double.PositiveInfinity, 1.0)] - public void ValidateInverseCumulativeDistribution(double x, double f) + public void ValidateInverseCumulativeDistribution(double x, double p) { var n = Normal.WithMeanStdDev(5.0, 2.0); - AssertHelpers.AlmostEqualRelative(x, n.InverseCumulativeDistribution(f), 14); - AssertHelpers.AlmostEqualRelative(x, Normal.InvCDF(5.0, 2.0, f), 14); + AssertHelpers.AlmostEqualRelative(x, n.InverseCumulativeDistribution(p), 14); + AssertHelpers.AlmostEqualRelative(x, Normal.InvCDF(5.0, 2.0, p), 14); } /// diff --git a/src/UnitTests/DistributionTests/Continuous/StudentTTests.cs b/src/UnitTests/DistributionTests/Continuous/StudentTTests.cs index 0f716cb6..8fb54c40 100644 --- a/src/UnitTests/DistributionTests/Continuous/StudentTTests.cs +++ b/src/UnitTests/DistributionTests/Continuous/StudentTTests.cs @@ -411,14 +411,6 @@ namespace MathNet.Numerics.UnitTests.DistributionTests.Continuous ied.Take(5).ToArray(); } - /// - /// Validate cumulative distribution. - /// - /// Location value. - /// Scale value. - /// Degrees of freedom. - /// Input X value. - /// Expected value. [TestCase(0.0, 1.0, 1.0, 0.0, 0.5)] [TestCase(0.0, 1.0, 1.0, 1.0, 0.75)] [TestCase(0.0, 1.0, 1.0, -1.0, 0.25)] @@ -432,10 +424,31 @@ namespace MathNet.Numerics.UnitTests.DistributionTests.Continuous [TestCase(0.0, 1.0, Double.PositiveInfinity, 0.0, 0.5)] [TestCase(0.0, 1.0, Double.PositiveInfinity, 1.0, 0.841344746068543)] [TestCase(0.0, 1.0, Double.PositiveInfinity, 2.0, 0.977249868051821)] - public void ValidateCumulativeDistribution(double location, double scale, double dof, double x, double c) + public void ValidateCumulativeDistribution(double location, double scale, double dof, double x, double p) { - var n = new StudentT(location, scale, dof); - AssertHelpers.AlmostEqualRelative(c, n.CumulativeDistribution(x), 13); + var dist = new StudentT(location, scale, dof); + Assert.That(dist.CumulativeDistribution(x), Is.EqualTo(p).Within(1e-13)); + Assert.That(StudentT.CDF(location, scale, dof, x), Is.EqualTo(p).Within(1e-13)); + } + + [TestCase(0.0, 1.0, 1.0, 0.0, 0.5)] + [TestCase(0.0, 1.0, 1.0, 1.0, 0.75)] + [TestCase(0.0, 1.0, 1.0, -1.0, 0.25)] + [TestCase(0.0, 1.0, 1.0, 2.0, 0.852416382349567)] + [TestCase(0.0, 1.0, 1.0, -2.0, 0.147583617650433)] + [TestCase(0.0, 1.0, 2.0, 0.0, 0.5)] + [TestCase(0.0, 1.0, 2.0, 1.0, 0.788675134594813)] + [TestCase(0.0, 1.0, 2.0, -1.0, 0.211324865405187)] + [TestCase(0.0, 1.0, 2.0, 2.0, 0.908248290463863)] + [TestCase(0.0, 1.0, 2.0, -2.0, 0.091751709536137)] + [TestCase(0.0, 1.0, Double.PositiveInfinity, 0.0, 0.5)] + [TestCase(0.0, 1.0, Double.PositiveInfinity, 1.0, 0.841344746068543)] + [TestCase(0.0, 1.0, Double.PositiveInfinity, 2.0, 0.977249868051821)] + public void ValidateInverseCumulativeDistribution(double location, double scale, double dof, double x, double p) + { + var dist = new StudentT(location, scale, dof); + Assert.That(dist.InverseCumulativeDistribution(p), Is.EqualTo(x).Within(1e-6)); + Assert.That(StudentT.InvCDF(location, scale, dof, p), Is.EqualTo(x).Within(1e-6)); } } }