diff --git a/src/Numerics/Distributions/Continuous/Normal.cs b/src/Numerics/Distributions/Continuous/Normal.cs
index d5363f6e..e8364fcd 100644
--- a/src/Numerics/Distributions/Continuous/Normal.cs
+++ b/src/Numerics/Distributions/Continuous/Normal.cs
@@ -282,6 +282,32 @@ namespace MathNet.Numerics.Distributions
get { return Double.PositiveInfinity; }
}
+ ///
+ /// Computes the density of the normal distribution.
+ ///
+ /// The mean of the normal distribution.
+ /// The standard deviation of the normal distribution.
+ /// The location at which to compute the density.
+ /// the density at .
+ 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);
+ }
+
+ ///
+ /// Computes the log density of the normal distribution.
+ ///
+ /// The mean of the normal distribution.
+ /// The standard deviation of the normal distribution.
+ /// The location at which to compute the density.
+ /// the log density at .
+ 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;
+ }
+
///
/// Computes the density of the normal distribution.
///
@@ -289,8 +315,7 @@ namespace MathNet.Numerics.Distributions
/// the density at .
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);
}
///
@@ -300,8 +325,7 @@ namespace MathNet.Numerics.Distributions
/// the log density at .
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);
}
///
diff --git a/src/Numerics/Distributions/Continuous/StudentT.cs b/src/Numerics/Distributions/Continuous/StudentT.cs
index b514c175..07d413b2 100644
--- a/src/Numerics/Distributions/Continuous/StudentT.cs
+++ b/src/Numerics/Distributions/Continuous/StudentT.cs
@@ -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
/// the density at .
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;
+ }
}
///
@@ -340,12 +355,19 @@ namespace MathNet.Numerics.Distributions
/// the log density at .
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);
+ }
}
///
@@ -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.
}
///
@@ -364,7 +387,7 @@ namespace MathNet.Numerics.Distributions
/// a sample from the distribution.
public double Sample()
{
- throw new NotImplementedException();
+ return _location + _scale * Sample(RandomSource, _dof);
}
///
@@ -373,7 +396,10 @@ namespace MathNet.Numerics.Distributions
/// a sequence of samples from the distribution.
public IEnumerable 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);
}
///
@@ -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);
+ }
+ }
+
+ ///
+ /// Samples standard student-t distributed random variables.
+ ///
+ /// The algorithm is method 2 in section 5, chapter 9
+ /// in L. Devroye's "Non-Uniform Random Variate Generation"
+ /// The random number generator to use.
+ /// The degrees of freedom for the standard student-t distribution.
+ /// a random number from the standard student-t distribution.
+ 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;
}
}
}
diff --git a/src/UnitTests/DistributionTests/CommonDistributionTests.cs b/src/UnitTests/DistributionTests/CommonDistributionTests.cs
index 7b4abc92..c3eebe91 100644
--- a/src/UnitTests/DistributionTests/CommonDistributionTests.cs
+++ b/src/UnitTests/DistributionTests/CommonDistributionTests.cs
@@ -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]
diff --git a/src/UnitTests/DistributionTests/Continuous/StudentTTests.cs b/src/UnitTests/DistributionTests/Continuous/StudentTTests.cs
index a5621fda..7b4dabb5 100644
--- a/src/UnitTests/DistributionTests/Continuous/StudentTTests.cs
+++ b/src/UnitTests/DistributionTests/Continuous/StudentTTests.cs
@@ -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(n.Mean, mean);
+ AssertEx.AreEqual(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(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(MathNet.Numerics.Constants.LogSqrt2PiE + Math.Log(n.StdDev), n.Entropy);
+ AssertEx.AreEqual(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(0.0, n.Skewness);
+ var n = new StudentT(location, scale, dof);
+ AssertEx.AreEqual(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(mean, n.Mode);
+ var n = new StudentT(location, scale, dof);
+ AssertEx.AreEqual(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(mean, n.Median);
+ var n = new StudentT(location, scale, dof);
+ AssertEx.AreEqual(location, n.Median);
}
[Test]
public void ValidateMinimum()
{
- var n = new Normal();
+ var n = new StudentT();
AssertEx.AreEqual(System.Double.NegativeInfinity, n.Minimum);
}
[Test]
public void ValidateMaximum()
{
- var n = new Normal();
+ var n = new StudentT();
AssertEx.AreEqual(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(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(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);
- }*/
}
}