From 666c02ff565acb62ec1c955cb09f1ffa02e21e1a Mon Sep 17 00:00:00 2001 From: Jurgen Van Gael Date: Wed, 5 May 2010 05:07:55 +0800 Subject: [PATCH] Fixed bug in NormalGamma distribution. Added documentation to StudentT distribution. --- src/Numerics/Distributions/Continuous/StudentT.cs | 14 ++++++++------ .../Distributions/Multivariate/NormalGamma.cs | 4 ++-- .../DistributionTests/Continuous/StudentTTests.cs | 8 ++++---- .../Multivariate/NormalGammaTests.cs | 5 +++-- 4 files changed, 17 insertions(+), 14 deletions(-) diff --git a/src/Numerics/Distributions/Continuous/StudentT.cs b/src/Numerics/Distributions/Continuous/StudentT.cs index a2b72107..ccdc6416 100644 --- a/src/Numerics/Distributions/Continuous/StudentT.cs +++ b/src/Numerics/Distributions/Continuous/StudentT.cs @@ -39,6 +39,8 @@ namespace MathNet.Numerics.Distributions /// We use a slightly generalized version (compared to Wikipedia) of the Student t-distribution. /// Namely, one which also parameterizes the location and scale. See the book "Bayesian Data Analysis" by Gelman /// et al. for more details. + /// The density of the Student t-distribution + /// p(x|mu,scale,dof) = Gamma((dof+1)/2) (1 + (x - mu)^2 / (scale * scale * dof))^(-(dof+1)/2) / (Gamma(dof/2)*Sqrt(dof*pi*scale)). /// The distribution will use the by default. /// Users can get/set the random number generator by using the property. /// The statistics classes will check all the incoming parameters whether they are in the allowed @@ -232,11 +234,11 @@ namespace MathNet.Numerics.Distributions { if (Double.IsPositiveInfinity(_dof)) { - return _scale; + return _scale * _scale; } else if (_dof > 2.0) { - return _dof * _scale / (_dof - 2.0); + return _dof * _scale * _scale / (_dof - 2.0); } else if (_dof > 1.0) { @@ -258,11 +260,11 @@ namespace MathNet.Numerics.Distributions { if (Double.IsPositiveInfinity(_dof)) { - return Math.Sqrt(_scale); + return Math.Sqrt(_scale * _scale); } else if (_dof > 2.0) { - return Math.Sqrt(_dof * _scale / (_dof - 2.0)); + return Math.Sqrt(_dof * _scale * _scale / (_dof - 2.0)); } else if (_dof > 1.0) { @@ -336,7 +338,7 @@ namespace MathNet.Numerics.Distributions // TODO JVG we can probably do a better job for Cauchy special case if (Double.IsPositiveInfinity(_dof)) { - return Normal.Density(_location, Math.Sqrt(_scale), x); + return Normal.Density(_location, _scale, x); } else { @@ -359,7 +361,7 @@ namespace MathNet.Numerics.Distributions // TODO JVG we can probably do a better job for Cauchy special case if (Double.IsPositiveInfinity(_dof)) { - return Normal.DensityLn(_location, Math.Sqrt(_scale), x); + return Normal.DensityLn(_location, _scale, x); } else { diff --git a/src/Numerics/Distributions/Multivariate/NormalGamma.cs b/src/Numerics/Distributions/Multivariate/NormalGamma.cs index 1da563ac..f3eb4485 100644 --- a/src/Numerics/Distributions/Multivariate/NormalGamma.cs +++ b/src/Numerics/Distributions/Multivariate/NormalGamma.cs @@ -242,11 +242,11 @@ namespace MathNet.Numerics.Distributions { if (Double.IsPositiveInfinity(_precisionInvScale)) { - return new StudentT(_meanLocation, _meanScale * _precisionShape, Double.PositiveInfinity); + return new StudentT(_meanLocation, 1.0 / (_meanScale * _precisionShape), Double.PositiveInfinity); } else { - return new StudentT(_meanLocation, _meanScale * _precisionShape / _precisionInvScale, 2.0 * _precisionShape); + return new StudentT(_meanLocation, Math.Sqrt(_precisionInvScale / (_meanScale * _precisionShape)), 2.0 * _precisionShape); } } diff --git a/src/UnitTests/DistributionTests/Continuous/StudentTTests.cs b/src/UnitTests/DistributionTests/Continuous/StudentTTests.cs index be234788..d4088451 100644 --- a/src/UnitTests/DistributionTests/Continuous/StudentTTests.cs +++ b/src/UnitTests/DistributionTests/Continuous/StudentTTests.cs @@ -172,8 +172,8 @@ namespace MathNet.Numerics.UnitTests.DistributionTests [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(0.0, 10.0, 2.5, 500.0)] + [Row(0.0, 10.0, Double.PositiveInfinity, 100.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)] @@ -190,8 +190,8 @@ namespace MathNet.Numerics.UnitTests.DistributionTests [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(0.0, 10.0, 2.5, 22.360679774997896964091736687313)] + [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, 2.2360679774997896964091736687313)] [Row(-5.0, 100.0, 1.5, Double.PositiveInfinity)] diff --git a/src/UnitTests/DistributionTests/Multivariate/NormalGammaTests.cs b/src/UnitTests/DistributionTests/Multivariate/NormalGammaTests.cs index 5777e610..aa244462 100644 --- a/src/UnitTests/DistributionTests/Multivariate/NormalGammaTests.cs +++ b/src/UnitTests/DistributionTests/Multivariate/NormalGammaTests.cs @@ -153,7 +153,7 @@ namespace MathNet.Numerics.UnitTests.DistributionTests [Test, MultipleAsserts] [Row(0.0, 1.0, 1.0, 1.0, 0.0, 1.0, 2.0)] [Row(10.0, 1.0, 2.0, 2.0, 10.0, 1.0, 4.0)] - [Row(10.0, 1.0, 2.0, Double.PositiveInfinity, 10.0, 2.0, Double.PositiveInfinity)] + [Row(10.0, 1.0, 2.0, Double.PositiveInfinity, 10.0, 0.5, Double.PositiveInfinity)] public void CanGetMeanMarginal(double meanLocation, double meanScale, double precShape, double precInvScale, double meanMarginalMean, double meanMarginalScale, double meanMarginalDoF) { @@ -209,7 +209,8 @@ namespace MathNet.Numerics.UnitTests.DistributionTests public void SampleFollowsCorrectDistribution() { Random rnd = new MersenneTwister(); - var cd = new NormalGamma(1.0, 4.0, 3.0, 3.5); + //var cd = new NormalGamma(1.0, 4.0, 3.0, 3.5); + var cd = new NormalGamma(1.0, 4.0, 7.0, 3.5); // Sample from the distribution. MeanPrecisionPair[] samples = new MeanPrecisionPair[CommonDistributionTests.NumberOfTestSamples];