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

Distributions: StudentT.InvCDF

optimization-3
Christoph Ruegg 13 years ago
parent
ccf4340ce0
commit
8f91ec929f
3 changed files with 75 additions and 19 deletions
Whitespace
Show all changes Ignore whitespace when comparing lines Ignore changes in amount of whitespace Ignore changes in whitespace at EOL
Split View
Diff Options
Show Stats Download Patch File Download Diff File
  1. 43
      src/Numerics/Distributions/StudentT.cs
  2. 16
      src/UnitTests/DistributionTests/Continuous/NormalTests.cs
  3. 35
      src/UnitTests/DistributionTests/Continuous/StudentTTests.cs

43
src/Numerics/Distributions/StudentT.cs
View File

@ -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);
}
/// <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"/>
/// <remarks>WARNING: currently not an explicit implementation, hence slow and unreliable.</remarks>
public double InverseCumulativeDistribution(double p)
{
return InvCDF(_location, _scale, _freedom, p);
}
/// <summary>
/// Samples student-t distributed random variables.
/// </summary>
@ -409,6 +423,35 @@ namespace MathNet.Numerics.Distributions
return x <= location ? ib : 1.0 - ib;
}
/// <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>
/// <param name="location">The location (μ) of the distribution.</param>
/// <param name="scale">The scale (σ) of the distribution. Range: σ > 0.</param>
/// <param name="freedom">The degrees of freedom (ν) for the distribution. Range: ν > 0.</param>
/// <returns>the inverse cumulative density at <paramref name="p"/>.</returns>
/// <seealso cref="InverseCumulativeDistribution"/>
/// <remarks>WARNING: currently not an explicit implementation, hence slow and unreliable.</remarks>
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);
}
/// <summary>
/// Generates a sample from the Student t-distribution.
/// </summary>

16
src/UnitTests/DistributionTests/Continuous/NormalTests.cs
View File

@ -458,7 +458,7 @@ namespace MathNet.Numerics.UnitTests.DistributionTests.Continuous
/// Validate cumulative distribution.
/// </summary>
/// <param name="x">Input X value.</param>
/// <param name="f">Expected value.</param>
/// <param name="p">Expected value.</param>
[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);
}
/// <summary>
/// Validate inverse cumulative distribution.
/// </summary>
/// <param name="x">Input X value.</param>
/// <param name="f">Expected value.</param>
/// <param name="p">Expected value.</param>
[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);
}
/// <summary>

35
src/UnitTests/DistributionTests/Continuous/StudentTTests.cs
View File

@ -411,14 +411,6 @@ namespace MathNet.Numerics.UnitTests.DistributionTests.Continuous
ied.Take(5).ToArray();
}
/// <summary>
/// Validate cumulative distribution.
/// </summary>
/// <param name="location">Location value.</param>
/// <param name="scale">Scale value.</param>
/// <param name="dof">Degrees of freedom.</param>
/// <param name="x">Input X value.</param>
/// <param name="c">Expected value.</param>
[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));
}
}
}

Write Preview
Loading…
Cancel
Save