diff --git a/src/Numerics/Distributions/Continuous/Beta.cs b/src/Numerics/Distributions/Continuous/Beta.cs index 31773bfe..bd0653a4 100644 --- a/src/Numerics/Distributions/Continuous/Beta.cs +++ b/src/Numerics/Distributions/Continuous/Beta.cs @@ -59,7 +59,7 @@ namespace MathNet.Numerics.Distributions private Random _random; /// - /// Initializes a new instance of the Beta distribution. + /// Initializes a new instance of the Beta class. /// /// The a shape parameter of the Beta distribution. /// The b shape parameter of the Beta distribution. @@ -73,6 +73,7 @@ namespace MathNet.Numerics.Distributions /// /// A string representation of the distribution. /// + /// A string representation of the Beta distribution. public override string ToString() { return "Beta(A = " + _shapeA + ", B = " + _shapeB + ")"; @@ -217,9 +218,9 @@ namespace MathNet.Numerics.Distributions get { return SpecialFunctions.BetaLn(_shapeA, _shapeB) - - (_shapeA - 1.0) * SpecialFunctions.DiGamma(_shapeA) - - (_shapeB - 1.0) * SpecialFunctions.DiGamma(_shapeB) - + (_shapeA + _shapeB - 2.0) * SpecialFunctions.DiGamma(_shapeA + _shapeB); + - ((_shapeA - 1.0) * SpecialFunctions.DiGamma(_shapeA)) + - ((_shapeB - 1.0) * SpecialFunctions.DiGamma(_shapeB)) + + ((_shapeA + _shapeB - 2.0) * SpecialFunctions.DiGamma(_shapeA + _shapeB)); } } @@ -505,8 +506,8 @@ namespace MathNet.Numerics.Distributions else { double a = SpecialFunctions.GammaLn(_shapeA + _shapeB) - SpecialFunctions.GammaLn(_shapeA) - SpecialFunctions.GammaLn(_shapeB); - double b = (x == 0.0 ? (_shapeA == 1.0 ? 0.0 : Double.NegativeInfinity) : (_shapeA - 1.0) * Math.Log(x)); - double c = (x == 1.0 ? (_shapeB == 1.0 ? 0.0 : Double.NegativeInfinity) : (_shapeB - 1.0) * Math.Log(1.0 - x)); + double b = x == 0.0 ? (_shapeA == 1.0 ? 0.0 : Double.NegativeInfinity) : (_shapeA - 1.0) * Math.Log(x); + double c = x == 1.0 ? (_shapeB == 1.0 ? 0.0 : Double.NegativeInfinity) : (_shapeB - 1.0) * Math.Log(1.0 - x); return a + b + c; } } diff --git a/src/Numerics/Distributions/Continuous/Gamma.cs b/src/Numerics/Distributions/Continuous/Gamma.cs index ea4cdb07..d5489377 100644 --- a/src/Numerics/Distributions/Continuous/Gamma.cs +++ b/src/Numerics/Distributions/Continuous/Gamma.cs @@ -67,7 +67,7 @@ namespace MathNet.Numerics.Distributions private Random _random; /// - /// Initializes a new instance of the Gamma distribution. + /// Initializes a new instance of the Gamma class. /// /// The shape of the Gamma distribution. /// The inverse scale of the Gamma distribution. @@ -304,7 +304,7 @@ namespace MathNet.Numerics.Distributions } else { - return _shape - Math.Log(_invScale) + SpecialFunctions.GammaLn(_shape) + (1.0 - _shape) * SpecialFunctions.DiGamma(_shape); + return _shape - Math.Log(_invScale) + SpecialFunctions.GammaLn(_shape) + ((1.0 - _shape) * SpecialFunctions.DiGamma(_shape)); } } } @@ -436,11 +436,11 @@ namespace MathNet.Numerics.Distributions } else if(_shape == 1.0) { - return Math.Log(_invScale) - _invScale*x; + return Math.Log(_invScale) - (_invScale * x); } else { - return _shape * Math.Log(_invScale) + (_shape - 1.0) * Math.Log(x) - _invScale * x - SpecialFunctions.GammaLn(_shape); + return (_shape * Math.Log(_invScale)) + ((_shape - 1.0) * Math.Log(x)) - (_invScale * x) - SpecialFunctions.GammaLn(_shape); } } @@ -554,25 +554,27 @@ namespace MathNet.Numerics.Distributions alphafix = System.Math.Pow(rnd.NextDouble(), 1.0 / shape); } - double d = a - 1.0 / 3.0; + double d = a - (1.0 / 3.0); double c = 1.0 / System.Math.Sqrt(9.0 * d); while (true) { double x = Normal.Sample(rnd, 0.0, 1.0); - double v = 1.0 + c * x; + double v = 1.0 + (c * x); while (v <= 0.0) { x = Normal.Sample(rnd, 0.0, 1.0); - v = 1.0 + c * x; + v = 1.0 + (c * x); } + v = v * v * v; double u = rnd.NextDouble(); x = x * x; - if (u < 1.0 - 0.0331 * x * x) + if (u < 1.0 - (0.0331 * x * x)) { return alphafix * d * v / invScale; } - if (System.Math.Log(u) < 0.5 * x + d * (1.0 - v + System.Math.Log(v))) + + if (System.Math.Log(u) < (0.5 * x) + (d * (1.0 - v + System.Math.Log(v)))) { return alphafix * d * v / invScale; } diff --git a/src/Numerics/SpecialFunctions.cs b/src/Numerics/SpecialFunctions.cs index 4761e779..5911581c 100644 --- a/src/Numerics/SpecialFunctions.cs +++ b/src/Numerics/SpecialFunctions.cs @@ -108,7 +108,8 @@ namespace MathNet.Numerics { s += Gamma_dk[i] / (i - z); } - return Constants.LnPi - Math.Log(Math.Sin(Math.PI * z)) - Math.Log(s) - Constants.LogTwoSqrtEOverPi - (0.5 - z) * Math.Log((0.5 - z + Gamma_r) / Math.E); + + return Constants.LnPi - Math.Log(Math.Sin(Math.PI * z)) - Math.Log(s) - Constants.LogTwoSqrtEOverPi - ((0.5 - z) * Math.Log((0.5 - z + Gamma_r) / Math.E)); } else { @@ -117,7 +118,8 @@ namespace MathNet.Numerics { s += Gamma_dk[i] / (z + i - 1.0); } - return Math.Log(s) + Constants.LogTwoSqrtEOverPi + (z - 0.5) * Math.Log((z - 0.5 + Gamma_r) / Math.E); + + return Math.Log(s) + Constants.LogTwoSqrtEOverPi + ((z - 0.5) * Math.Log((z - 0.5 + Gamma_r) / Math.E)); } } @@ -144,6 +146,7 @@ namespace MathNet.Numerics { s += Gamma_dk[i] / (i - z); } + return Math.PI / (Math.Sin(Math.PI * z) * s * Constants.TwoSqrtEOverPi * Math.Pow((0.5 - z + Gamma_r) / Math.E, 0.5 - z)); } else @@ -153,6 +156,7 @@ namespace MathNet.Numerics { s += Gamma_dk[i] / (z + i - 1.0); } + return s * Constants.TwoSqrtEOverPi * Math.Pow((z - 0.5 + Gamma_r) / Math.E, z - 0.5); } } @@ -195,12 +199,12 @@ namespace MathNet.Numerics // Use inversion formula for negative numbers. if (x < 0) { - return DiGamma(1.0 - x) + System.Math.PI / System.Math.Tan(-System.Math.PI * x); + return DiGamma(1.0 - x) + (System.Math.PI / System.Math.Tan(-System.Math.PI * x)); } if (x <= s) { - return d1 - 1 / x + d2 * x; + return d1 - (1 / x) + (d2 * x); } double result = 0; @@ -213,10 +217,10 @@ namespace MathNet.Numerics if (x >= c) { double r = 1 / x; - result += System.Math.Log(x) - 0.5 * r; + result += System.Math.Log(x) - (0.5 * r); r *= r; - result -= r * (s3 - r * (s4 - r * (s5 - r * (s6 - r * s7)))); + result -= r * (s3 - (r * (s4 - (r * (s5 - (r * (s6 - (r * s7)))))))); } return result; @@ -226,6 +230,7 @@ namespace MathNet.Numerics { throw new NotImplementedException(); } + public static double BetaLn(double a, double b) { throw new NotImplementedException();