From f8792b77f40f5f53b46e3bcfde319c215974e127 Mon Sep 17 00:00:00 2001 From: Christoph Ruegg Date: Wed, 29 Aug 2012 11:56:20 +0200 Subject: [PATCH] Dist: StudentT density more robust for very large degrees of freedom #44 --- src/Numerics/Distributions/Continuous/StudentT.cs | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/src/Numerics/Distributions/Continuous/StudentT.cs b/src/Numerics/Distributions/Continuous/StudentT.cs index 072050a1..53113176 100644 --- a/src/Numerics/Distributions/Continuous/StudentT.cs +++ b/src/Numerics/Distributions/Continuous/StudentT.cs @@ -352,15 +352,14 @@ namespace MathNet.Numerics.Distributions public double Density(double x) { // TODO JVG we can probably do a better job for Cauchy special case - if (Double.IsPositiveInfinity(_dof)) + if (_dof >= 1e+8d) { return Normal.Density(_location, _scale, x); } var d = (x - _location) / _scale; - return SpecialFunctions.Gamma((_dof + 1.0) / 2.0) + return Math.Exp(SpecialFunctions.GammaLn((_dof + 1.0) / 2.0) - SpecialFunctions.GammaLn(_dof / 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; } @@ -373,7 +372,7 @@ namespace MathNet.Numerics.Distributions public double DensityLn(double x) { // TODO JVG we can probably do a better job for Cauchy special case - if (Double.IsPositiveInfinity(_dof)) + if (_dof >= 1e+8d) { return Normal.DensityLn(_location, _scale, x); }