Distributions: Chi, ChiSquared and NormalGamma more robust for large shape parameter

src/Numerics/Distributions/Chi.cs

@@ -317,6 +317,11 @@ namespace MathNet.Numerics.Distributions
                 throw new ArgumentException(Resources.InvalidDistributionParameters);
             }
 
+            if (freedom > 160.0)
+            {
+                return Math.Exp(PDFLn(freedom, x));
+            }
+
             return (Math.Pow(2.0, 1.0 - (freedom/2.0))*Math.Pow(x, freedom - 1.0)*Math.Exp(-x*x/2.0))/SpecialFunctions.Gamma(freedom/2.0);
         }
 
@@ -351,7 +356,7 @@ namespace MathNet.Numerics.Distributions
                 throw new ArgumentException(Resources.InvalidDistributionParameters);
             }
 
-            return SpecialFunctions.GammaLowerIncomplete(freedom/2.0, x*x/2.0)/SpecialFunctions.Gamma(freedom/2.0);
+            return SpecialFunctions.GammaLowerRegularized(freedom/2.0, x*x/2.0);
         }
 
         /// <summary>

src/Numerics/Distributions/ChiSquared.cs

@@ -322,6 +322,11 @@ namespace MathNet.Numerics.Distributions
                 throw new ArgumentException(Resources.InvalidDistributionParameters);
             }
 
+            if (freedom > 160.0)
+            {
+                return Math.Exp(PDFLn(freedom, x));
+            }
+
             return (Math.Pow(x, (freedom/2.0) - 1.0)*Math.Exp(-x/2.0))/(Math.Pow(2.0, freedom/2.0)*SpecialFunctions.Gamma(freedom/2.0));
         }
 

src/Numerics/Distributions/NormalGamma.cs

@@ -283,6 +283,11 @@ namespace MathNet.Numerics.Distributions
                 throw new NotSupportedException();
             }
 
+            if (_precisionShape > 160.0)
+            {
+                return Math.Exp(DensityLn(mean, prec));
+            }
+
             // double e = -0.5 * prec * (mean - _meanLocation) * (mean - _meanLocation) - prec * _precisionInvScale;
             // return Math.Pow(prec * _precisionInvScale, _precisionShape) * Math.Exp(e) / (Constants.Sqrt2Pi * Math.Sqrt(prec) * SpecialFunctions.Gamma(_precisionShape));
             double e = -(0.5*prec*_meanScale*(mean - _meanLocation)*(mean - _meanLocation)) - (prec*_precisionInvScale);