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

Added more unit tests for the StudentT distribution. 

Added sampling methods for the StudentT distribution.
la-knuth
Jurgen Van Gael 17 years ago
parent
c66021cd15
commit
e3d4c58cdc
4 changed files with 207 additions and 160 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. 32
      src/Numerics/Distributions/Continuous/Normal.cs
  2. 90
      src/Numerics/Distributions/Continuous/StudentT.cs
  3. 1
      src/UnitTests/DistributionTests/CommonDistributionTests.cs
  4. 244
      src/UnitTests/DistributionTests/Continuous/StudentTTests.cs

32
src/Numerics/Distributions/Continuous/Normal.cs
View File

@ -282,6 +282,32 @@ namespace MathNet.Numerics.Distributions
get { return Double.PositiveInfinity; }
}
/// <summary>
/// Computes the density of the normal distribution.
/// </summary>
/// <param name="mean">The mean of the normal distribution.</param>
/// <param name="sdev">The standard deviation of the normal distribution.</param>
/// <param name="x">The location at which to compute the density.</param>
/// <returns>the density at <paramref name="x"/>.</returns>
internal static double Density(double mean, double sdev, double x)
{
double d = (x - mean) / sdev;
return Math.Exp(-0.5 * d * d) / (Constants.Sqrt2Pi * sdev);
}
/// <summary>
/// Computes the log density of the normal distribution.
/// </summary>
/// <param name="mean">The mean of the normal distribution.</param>
/// <param name="sdev">The standard deviation of the normal distribution.</param>
/// <param name="x">The location at which to compute the density.</param>
/// <returns>the log density at <paramref name="x"/>.</returns>
internal static double DensityLn(double mean, double sdev, double x)
{
double d = (x - mean) / sdev;
return (-0.5 * d * d) - Math.Log(sdev) - Constants.LogSqrt2Pi;
}
/// <summary>
/// Computes the density of the normal distribution.
/// </summary>
@ -289,8 +315,7 @@ namespace MathNet.Numerics.Distributions
/// <returns>the density at <paramref name="x"/>.</returns>
public double Density(double x)
{
double d = (x - _mean) / _stdDev;
return Math.Exp(-0.5 * d * d) / (Constants.Sqrt2Pi * _stdDev);
return Density(_mean, _stdDev, x);
}
/// <summary>
@ -300,8 +325,7 @@ namespace MathNet.Numerics.Distributions
/// <returns>the log density at <paramref name="x"/>.</returns>
public double DensityLn(double x)
{
double d = (x - _mean) / _stdDev;
return (-0.5 * d * d) - Math.Log(_stdDev) - Constants.LogSqrt2Pi;
return DensityLn(_mean, _stdDev, x);
}
/// <summary>

90
src/Numerics/Distributions/Continuous/StudentT.cs
View File

@ -230,9 +230,13 @@ namespace MathNet.Numerics.Distributions
{
get
{
if (_dof > 2.0)
if (Double.IsPositiveInfinity(_dof))
{
return _dof / (_dof - 2.0) / _scale;
return _scale;
}
else if (_dof > 2.0)
{
return _dof * _scale / (_dof - 2.0);
}
else if (_dof > 1.0)
{
@ -240,7 +244,7 @@ namespace MathNet.Numerics.Distributions
}
else
{
throw new Exception(Resources.UndefinedMoment);
return Double.NaN;
}
}
}
@ -252,9 +256,13 @@ namespace MathNet.Numerics.Distributions
{
get
{
if (_dof > 2.0)
if (Double.IsPositiveInfinity(_dof))
{
return Math.Sqrt(_dof / (_dof - 2.0));
return Math.Sqrt(_scale);
}
else if (_dof > 2.0)
{
return Math.Sqrt(_dof * _scale / (_dof - 2.0));
}
else if (_dof > 1.0)
{
@ -262,7 +270,7 @@ namespace MathNet.Numerics.Distributions
}
else
{
throw new Exception(Resources.UndefinedMoment);
return Double.NaN;
}
}
}
@ -325,12 +333,19 @@ namespace MathNet.Numerics.Distributions
/// <returns>the density at <paramref name="x"/>.</returns>
public double Density(double x)
{
double d = (x - _location) / _scale;
return SpecialFunctions.Gamma((_dof + 1.0) / 2.0)
* Math.Pow(1.0 + d * d / _dof, -0.5 * (_dof + 1.0))
/ SpecialFunctions.Gamma(_dof / 2.0)
/ Math.Sqrt(_dof * Math.PI)
/ _scale;
if (Double.IsPositiveInfinity(_dof))
{
return Normal.Density(_location, Math.Sqrt(_scale), x);
}
else
{
double d = (x - _location) / _scale;
return SpecialFunctions.Gamma((_dof + 1.0) / 2.0)
* Math.Pow(1.0 + d * d / _dof, -0.5 * (_dof + 1.0))
/ SpecialFunctions.Gamma(_dof / 2.0)
/ Math.Sqrt(_dof * Math.PI)
/ _scale;
}
}
/// <summary>
@ -340,12 +355,19 @@ namespace MathNet.Numerics.Distributions
/// <returns>the log density at <paramref name="x"/>.</returns>
public double DensityLn(double x)
{
double d = (x - _location) / _scale;
return SpecialFunctions.GammaLn((_dof + 1.0) / 2.0)
- 0.5 * (_dof + 1.0) * Math.Log(1.0 + d * d / _dof)
- SpecialFunctions.GammaLn(_dof / 2.0)
-0.5 * Math.Log(_dof * Math.PI)
- Math.Log(_scale);
if (Double.IsPositiveInfinity(_dof))
{
return Normal.DensityLn(_location, Math.Sqrt(_scale), x);
}
else
{
double d = (x - _location) / _scale;
return SpecialFunctions.GammaLn((_dof + 1.0) / 2.0)
- 0.5 * (_dof + 1.0) * Math.Log(1.0 + d * d / _dof)
- SpecialFunctions.GammaLn(_dof / 2.0)
- 0.5 * Math.Log(_dof * Math.PI)
- Math.Log(_scale);
}
}
/// <summary>
@ -356,6 +378,7 @@ namespace MathNet.Numerics.Distributions
public double CumulativeDistribution(double x)
{
throw new NotImplementedException();
// TODO Jurgen: once this is implemented; enable the StudentT stuff in commondistributiontests.
}
/// <summary>
@ -364,7 +387,7 @@ namespace MathNet.Numerics.Distributions
/// <returns>a sample from the distribution.</returns>
public double Sample()
{
throw new NotImplementedException();
return _location + _scale * Sample(RandomSource, _dof);
}
/// <summary>
@ -373,7 +396,10 @@ namespace MathNet.Numerics.Distributions
/// <returns>a sequence of samples from the distribution.</returns>
public IEnumerable<double> Samples()
{
throw new NotImplementedException();
while (true)
{
yield return _location + _scale * Sample(RandomSource, _dof);
}
}
#endregion
@ -392,7 +418,7 @@ namespace MathNet.Numerics.Distributions
throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters);
}
throw new NotImplementedException();
return location + scale * Sample(rng, dof);
}
/// <summary>
@ -410,7 +436,27 @@ namespace MathNet.Numerics.Distributions
throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters);
}
throw new NotImplementedException();
while (true)
{
yield return location + scale * Sample(rng, dof);
}
}
/// <summary>
/// Samples standard student-t distributed random variables.
/// </summary>
/// <remarks>The algorithm is method 2 in section 5, chapter 9
/// in L. Devroye's "Non-Uniform Random Variate Generation"</remarks>
/// <param name="rnd">The random number generator to use.</param>
/// <param name="dof">The degrees of freedom for the standard student-t distribution.</param>
/// <returns>a random number from the standard student-t distribution.</returns>
internal static double Sample(Random rnd, double dof)
{
double dummy = 0.0;
var n = Normal.SampleBoxMuller(rnd, out dummy);
var g = Gamma.Sample(rnd, dof / 2.0, 1.0);
return Math.Sqrt(2.0 * dof / g) * n;
}
}
}

1
src/UnitTests/DistributionTests/CommonDistributionTests.cs
View File

@ -66,6 +66,7 @@ namespace MathNet.Numerics.UnitTests.DistributionTests
continuousDistributions.Add(new Normal(0.0, 1.0));
continuousDistributions.Add(new Weibull(1.0, 1.0));
continuousDistributions.Add(new LogNormal(1.0, 1.0));
//continuousDistributions.Add(new StudentT(0.0, 1.0, 3.0));
}
[Test]

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

@ -159,214 +159,190 @@ namespace MathNet.Numerics.UnitTests.DistributionTests
[Row(0.0, 10.0, Double.PositiveInfinity, 0.0)]
[Row(10.0, 1.0, 1.0, Double.NaN)]
[Row(-5.0, 100.0, 1.5, -5.0)]
[Row(0.0, Double.PositiveInfinity, 1.0)]
[Row(0.0, Double.PositiveInfinity, 1.0, Double.NaN)]
public void ValidateMean(double location, double scale, double dof, double mean)
{
var n = new StudentT(location, scale, dof);
AssertEx.AreEqual<double>(n.Mean, mean);
AssertEx.AreEqual<double>(mean, n.Mean);
}
/*
[Test]
[Row(0.0, 1.0, 1.0)]
[Row(0.0, 0.1, 1.0)]
[Row(0.0, 1.0, 3.0)]
[Row(0.0, 10.0, 1.0)]
[Row(0.0, 10.0, 2.0)]
[Row(0.0, 10.0, 3.0)]
[Row(0.0, 10.0, Double.PositiveInfinity)]
[Row(10.0, 1.0, 1.0)]
[Row(-5.0, 100.0, 1.0)]
[Row(0.0, Double.PositiveInfinity, 1.0)]
[Row(0.0, 1.0, 1.0, Double.NaN)]
[Row(0.0, 0.1, 1.0, Double.NaN)]
[Row(0.0, 1.0, 3.0, 3.0)]
[Row(0.0, 10.0, 1.0, Double.NaN)]
[Row(0.0, 10.0, 2.0, Double.PositiveInfinity)]
[Row(0.0, 10.0, 2.5, 50.0)]
[Row(0.0, 10.0, Double.PositiveInfinity, 10.0)]
[Row(10.0, 1.0, 1.0, Double.NaN)]
[Row(10.0, 1.0, 2.5, 5.0)]
[Row(-5.0, 100.0, 1.5, Double.PositiveInfinity)]
[Row(0.0, Double.PositiveInfinity, 1.0, Double.NaN)]
public void ValidateVariance(double location, double scale, double dof, double var)
{
var n = new StudentT(location, scale, dof);
AssertEx.AreEqual<double>(n.Variance, location);
}
[Test]
[Row(-0.0)]
[Row(0.0)]
[Row(0.1)]
[Row(1.0)]
[Row(10.0)]
[Row(Double.PositiveInfinity)]
public void Entropy(double sdev)
{
var n = new Normal(1.0, sdev);
AssertEx.AreEqual<double>(MathNet.Numerics.Constants.LogSqrt2PiE + Math.Log(n.StdDev), n.Entropy);
AssertEx.AreEqual<double>(var, n.Variance);
}
[Test]
[Row(-0.0)]
[Row(0.0)]
[Row(0.1)]
[Row(1.0)]
[Row(10.0)]
[Row(Double.PositiveInfinity)]
public void ValidateSkewness(double sdev)
[Row(0.0, 1.0, 1.0, Double.NaN)]
[Row(0.0, 0.1, 1.0, Double.NaN)]
[Row(0.0, 1.0, 3.0, 1.7320508075688772935274463415059)]
[Row(0.0, 10.0, 1.0, Double.NaN)]
[Row(0.0, 10.0, 2.0, Double.PositiveInfinity)]
[Row(0.0, 10.0, 2.5, 7.0710678118654752440084436210485)]
[Row(0.0, 10.0, Double.PositiveInfinity, 3.1622776601683793319988935444327)]
[Row(10.0, 1.0, 1.0, Double.NaN)]
[Row(10.0, 1.0, 2.5, 2.2360679774997896964091736687313)]
[Row(-5.0, 100.0, 1.5, Double.PositiveInfinity)]
[Row(0.0, Double.PositiveInfinity, 1.0, Double.NaN)]
public void ValidateStdDev(double location, double scale, double dof, double sdev)
{
var n = new Normal(1.0, sdev);
AssertEx.AreEqual<double>(0.0, n.Skewness);
var n = new StudentT(location, scale, dof);
AssertEx.AreEqual<double>(sdev, n.StdDev);
}
[Test]
[Row(Double.NegativeInfinity)]
[Row(-0.0)]
[Row(0.0)]
[Row(0.1)]
[Row(1.0)]
[Row(10.0)]
[Row(Double.PositiveInfinity)]
public void ValidateMode(double mean)
[Row(0.0, 1.0, 1.0)]
[Row(0.0, 0.1, 1.0)]
[Row(0.0, 1.0, 3.0)]
[Row(0.0, 10.0, 1.0)]
[Row(0.0, 10.0, 2.0)]
[Row(0.0, 10.0, 2.5)]
[Row(0.0, 10.0, Double.PositiveInfinity)]
[Row(10.0, 1.0, 1.0)]
[Row(10.0, 1.0, 2.5)]
[Row(-5.0, 100.0, 1.5)]
[Row(0.0, Double.PositiveInfinity, 1.0)]
public void ValidateMode(double location, double scale, double dof)
{
var n = new Normal(mean, 1.0);
AssertEx.AreEqual<double>(mean, n.Mode);
var n = new StudentT(location, scale, dof);
AssertEx.AreEqual<double>(location, n.Mode);
}
[Test]
[Row(Double.NegativeInfinity)]
[Row(-0.0)]
[Row(0.0)]
[Row(0.1)]
[Row(1.0)]
[Row(10.0)]
[Row(Double.PositiveInfinity)]
public void ValidateMedian(double mean)
[Row(0.0, 1.0, 1.0)]
[Row(0.0, 0.1, 1.0)]
[Row(0.0, 1.0, 3.0)]
[Row(0.0, 10.0, 1.0)]
[Row(0.0, 10.0, 2.0)]
[Row(0.0, 10.0, 2.5)]
[Row(0.0, 10.0, Double.PositiveInfinity)]
[Row(10.0, 1.0, 1.0)]
[Row(10.0, 1.0, 2.5)]
[Row(-5.0, 100.0, 1.5)]
[Row(0.0, Double.PositiveInfinity, 1.0)]
public void ValidateMedian(double location, double scale, double dof)
{
var n = new Normal(mean, 1.0);
AssertEx.AreEqual<double>(mean, n.Median);
var n = new StudentT(location, scale, dof);
AssertEx.AreEqual<double>(location, n.Median);
}
[Test]
public void ValidateMinimum()
{
var n = new Normal();
var n = new StudentT();
AssertEx.AreEqual<double>(System.Double.NegativeInfinity, n.Minimum);
}
[Test]
public void ValidateMaximum()
{
var n = new Normal();
var n = new StudentT();
AssertEx.AreEqual<double>(System.Double.PositiveInfinity, n.Maximum);
}
[Test]
[Row(0.0, 0.0)]
[Row(0.0, 0.1)]
[Row(0.0, 1.0)]
[Row(0.0, 10.0)]
[Row(10.0, 1.0)]
[Row(-5.0, 100.0)]
[Row(0.0, Double.PositiveInfinity)]
public void ValidateDensity(double mean, double sdev)
[Row(0.0, 1.0, 1.0, 0.0, 0.318309886183791)]
[Row(0.0, 1.0, 1.0, 1.0, 0.159154943091895)]
[Row(0.0, 1.0, 1.0, -1.0, 0.159154943091895)]
[Row(0.0, 1.0, 1.0, 2.0, 0.063661977236758)]
[Row(0.0, 1.0, 1.0, -2.0, 0.063661977236758)]
[Row(0.0, 1.0, 2.0, 0.0, 0.353553390593274)]
[Row(0.0, 1.0, 2.0, 1.0, 0.192450089729875)]
[Row(0.0, 1.0, 2.0, -1.0, 0.192450089729875)]
[Row(0.0, 1.0, 2.0, 2.0, 0.068041381743977)]
[Row(0.0, 1.0, 2.0, -2.0, 0.068041381743977)]
[Row(0.0, 1.0, Double.PositiveInfinity, 0.0, 0.398942280401433)]
[Row(0.0, 1.0, Double.PositiveInfinity, 1.0, 0.241970724519143)]
[Row(0.0, 1.0, Double.PositiveInfinity, 2.0, 0.053990966513188)]
public void ValidateDensity(double location, double scale, double dof, double x, double p)
{
var n = Normal.WithMeanStdDev(mean, sdev);
for(int i = 0; i < 11; i++)
{
double x = i - 5.0;
double d = (mean - x)/sdev;
double pdf = Math.Exp(-0.5*d*d)/(sdev*Constants.Sqrt2Pi);
AssertEx.AreEqual<double>(pdf, n.Density(x));
}
var n = new StudentT(location, scale, dof);
AssertHelpers.AlmostEqual(p, n.Density(x), 13);
}
[Test]
[Row(0.0, 0.0)]
[Row(0.0, 0.1)]
[Row(0.0, 1.0)]
[Row(0.0, 10.0)]
[Row(10.0, 1.0)]
[Row(-5.0, 100.0)]
[Row(0.0, Double.PositiveInfinity)]
public void ValidateDensityLn(double mean, double sdev)
[Row(0.0, 1.0, 1.0, 0.0, -1.144729885849399)]
[Row(0.0, 1.0, 1.0, 1.0, -1.837877066409348)]
[Row(0.0, 1.0, 1.0, -1.0, -1.837877066409348)]
[Row(0.0, 1.0, 1.0, 2.0, -2.754167798283503)]
[Row(0.0, 1.0, 1.0, -2.0, -2.754167798283503)]
[Row(0.0, 1.0, 2.0, 0.0, -1.039720770839917)]
[Row(0.0, 1.0, 2.0, 1.0, -1.647918433002166)]
[Row(0.0, 1.0, 2.0, -1.0, -1.647918433002166)]
[Row(0.0, 1.0, 2.0, 2.0, -2.687639203842085)]
[Row(0.0, 1.0, 2.0, -2.0, -2.687639203842085)]
[Row(0.0, 1.0, Double.PositiveInfinity, 0.0, -0.918938533204672)]
[Row(0.0, 1.0, Double.PositiveInfinity, 1.0, -1.418938533204674)]
[Row(0.0, 1.0, Double.PositiveInfinity, 2.0, -2.918938533204674)]
public void ValidateDensityLn(double location, double scale, double dof, double x, double p)
{
var n = Normal.WithMeanStdDev(mean, sdev);
for (int i = 0; i < 11; i++)
{
double x = i - 5.0;
double d = (mean - x) / sdev;
double pdfln = -0.5 * d * d - Math.Log(sdev) - Constants.LogSqrt2Pi;
AssertEx.AreEqual<double>(pdfln, n.DensityLn(x));
}
var n = new StudentT(location, scale, dof);
AssertHelpers.AlmostEqual(p, n.DensityLn(x), 13);
}
[Test]
public void CanSampleStatic()
{
var d = Normal.Sample(new Random(), 0.0, 1.0);
var d = StudentT.Sample(new Random(), 0.0, 1.0, 3.0);
}
[Test]
public void CanSampleSequenceStatic()
{
var ied = Normal.Samples(new Random(), 0.0, 1.0);
var ied = StudentT.Samples(new Random(), 0.0, 1.0, 3.0);
var arr = ied.Take(5).ToArray();
}
[Test]
[ExpectedException(typeof(ArgumentOutOfRangeException))]
public void FailSampleStatic()
[Row(0.0, Double.NaN, 1.0)]
[Row(0.0, 1.0, Double.NaN)]
[Row(0.0, -1.0, 1.0)]
[Row(0.0, 1.0, -1.0)]
[Row(Double.NaN, 1.0, Double.NaN)]
public void FailSampleStatic(double location, double scale, double dof)
{
var d = Normal.Sample(new Random(), 0.0, -1.0);
var d = StudentT.Sample(new Random(), location, scale, dof);
}
[Test]
[ExpectedException(typeof(ArgumentOutOfRangeException))]
public void FailSampleSequenceStatic()
[Row(0.0, Double.NaN, 1.0)]
[Row(0.0, 1.0, Double.NaN)]
[Row(0.0, -1.0, 1.0)]
[Row(0.0, 1.0, -1.0)]
[Row(Double.NaN, 1.0, 1.0)]
public void FailSampleSequenceStatic(double location, double scale, double dof)
{
var ied = Normal.Samples(new Random(), 0.0, -1.0).First();
var ied = StudentT.Samples(new Random(), location, scale, dof);
var e = ied.Take(5).ToArray();
}
[Test]
public void CanSample()
{
var n = new Normal();
var n = new StudentT();
var d = n.Sample();
}
[Test]
public void CanSampleSequence()
{
var n = new Normal();
var n = new StudentT();
var ied = n.Samples();
var e = ied.Take(5).ToArray();
}
[Test]
[Row(Double.NegativeInfinity, 0.0)]
[Row(-5.0, 0.00000028665157187919391167375233287464535385442301361187883)]
[Row(-2.0, 0.0002326290790355250363499258867279847735487493358890356)]
[Row(-0.0, 0.0062096653257761351669781045741922211278977469230927036)]
[Row(0.0, 0.0062096653257761351669781045741922211278977469230927036)]
[Row(4.0, 0.30853753872598689636229538939166226011639782444542207)]
[Row(5.0, 0.5)]
[Row(6.0, 0.69146246127401310363770461060833773988360217555457859)]
[Row(10.0, 0.9937903346742238648330218954258077788721022530769078)]
[Row(Double.PositiveInfinity, 1.0)]
public void ValidateCumulativeDistribution(double x, double f)
{
var n = Normal.WithMeanStdDev(5.0, 2.0);
AssertHelpers.AlmostEqual(f, n.CumulativeDistribution(x), 10);
}
[Test]
[Row(Double.NegativeInfinity, 0.0)]
[Row(-5.0, 0.00000028665157187919391167375233287464535385442301361187883)]
[Row(-2.0, 0.0002326290790355250363499258867279847735487493358890356)]
[Row(-0.0, 0.0062096653257761351669781045741922211278977469230927036)]
[Row(0.0, 0.0062096653257761351669781045741922211278977469230927036)]
[Row(4.0, 0.30853753872598689636229538939166226011639782444542207)]
[Row(5.0, 0.5)]
[Row(6.0, 0.69146246127401310363770461060833773988360217555457859)]
[Row(10.0, 0.9937903346742238648330218954258077788721022530769078)]
[Row(Double.PositiveInfinity, 1.0)]
public void ValidateInverseCumulativeDistribution(double x, double f)
{
var n = Normal.WithMeanStdDev(5.0, 2.0);
AssertHelpers.AlmostEqual(x, n.InverseCumulativeDistribution(f), 15);
}*/
}
}

Write Preview
Loading…
Cancel
Save