diff --git a/src/Managed.UnitTests/AssertHelpers.cs b/src/Managed.UnitTests/AssertHelpers.cs index bef96af9..5d3c397e 100644 --- a/src/Managed.UnitTests/AssertHelpers.cs +++ b/src/Managed.UnitTests/AssertHelpers.cs @@ -38,13 +38,19 @@ namespace MathNet.Numerics.UnitTests class AssertHelpers { /// - /// Asserts that the expected value and the actual value are equal up to a certain number of decimal places. + /// Asserts that the expected value and the actual value are equal up to a certain number of decimal places. If both + /// and are NaN then no assert is thrown. /// /// The expected value. /// The actual value. /// The number of decimal places to agree on. public static void AlmostEqual(double expected, double actual, int decimalPlaces) { + if(double.IsNaN(expected) && double.IsNaN(actual)) + { + return; + } + bool pass = Precision.AlmostEqualInDecimalPlaces(expected, actual, decimalPlaces); if (!pass) { diff --git a/src/Managed.UnitTests/DistributionTests/Continuous/BetaTests.cs b/src/Managed.UnitTests/DistributionTests/Continuous/BetaTests.cs index ebc129eb..abab4984 100644 --- a/src/Managed.UnitTests/DistributionTests/Continuous/BetaTests.cs +++ b/src/Managed.UnitTests/DistributionTests/Continuous/BetaTests.cs @@ -124,22 +124,22 @@ namespace MathNet.Numerics.UnitTests.DistributionTests [Test] [Row(0.0, 0.0, 0.5)] - [Row(0.0, 0.1, 0.1)] + [Row(0.0, 0.1, 0.0)] [Row(1.0, 0.0, 1.0)] [Row(1.0, 1.0, 0.5)] [Row(9.0, 1.0, 0.9)] [Row(5.0, 100.0, 0.047619047619047619047616)] - [Row(1.0, Double.PositiveInfinity, 1.0)] - [Row(Double.PositiveInfinity, 1.0, 0.0)] - [Row(0.0, Double.PositiveInfinity, 1.0)] - [Row(Double.PositiveInfinity, 0.0, 0.0)] + [Row(1.0, Double.PositiveInfinity, 0.0)] + [Row(Double.PositiveInfinity, 1.0, 1.0)] + [Row(0.0, Double.PositiveInfinity, 0.0)] + [Row(Double.PositiveInfinity, 0.0, 1.0)] public void ValidateMean(double a, double b, double mean) { var n = new Beta(a, b); AssertEx.AreEqual(mean, n.Mean); } - [Test] + [Test, Ignore("Depending on Special Functions")] [Row(0.0, 0.0, 0.5)] [Row(0.0, 0.1, 0.1)] [Row(1.0, 0.0, 1.0)] @@ -175,8 +175,8 @@ namespace MathNet.Numerics.UnitTests.DistributionTests [Test] [Row(0.0, 0.0, 0.5)] - [Row(0.0, 0.1, 1.0)] - [Row(1.0, 0.0, 0.0)] + [Row(0.0, 0.1, 0.0)] + [Row(1.0, 0.0, 1.0)] [Row(1.0, 1.0, 0.5)] [Row(9.0, 1.0, 1.0)] [Row(5.0, 100.0, 0.038834951456310676243255386452801758423447608947753906)] @@ -338,7 +338,7 @@ namespace MathNet.Numerics.UnitTests.DistributionTests AssertHelpers.AlmostEqual(pdfln, n.DensityLn(x), 15); } - [Test] + [Test, Ignore("Depending on Special Functions")] [Row(0.0, 0.0, 0.0, 0.5)] [Row(0.0, 0.0, 0.5, 0.5)] [Row(0.0, 0.0, 1.0, 1.0)] diff --git a/src/Managed.UnitTests/DistributionTests/Continuous/GammaTests.cs b/src/Managed.UnitTests/DistributionTests/Continuous/GammaTests.cs index ae7557d4..e4d6b1a4 100644 --- a/src/Managed.UnitTests/DistributionTests/Continuous/GammaTests.cs +++ b/src/Managed.UnitTests/DistributionTests/Continuous/GammaTests.cs @@ -168,13 +168,13 @@ namespace MathNet.Numerics.UnitTests.DistributionTests } [Test] - [Row(0.0, 0.0, 0.0)] + [Row(0.0, 0.0, Double.NaN)] [Row(1.0, 0.1, 10.0)] [Row(1.0, 1.0, 1.0)] [Row(10.0, 10.0, 1.0)] [Row(10.0, 1.0, 10.0)] - [Row(10.0, Double.PositiveInfinity, 0.0)] - public void CanGetMean(double shape, double invScale, double mean) + [Row(10.0, Double.PositiveInfinity, 10.0)] + public void ValidateMean(double shape, double invScale, double mean) { var n = new Gamma(shape, invScale); AssertEx.AreEqual(mean, n.Mean); @@ -187,7 +187,7 @@ namespace MathNet.Numerics.UnitTests.DistributionTests [Row(10.0, 10.0, 0.1)] [Row(10.0, 1.0, 10.0)] [Row(10.0, Double.PositiveInfinity, 0.0)] - public void CanGetVariance(double shape, double invScale, double var) + public void ValidateVariance(double shape, double invScale, double var) { var n = new Gamma(shape, invScale); AssertEx.AreEqual(var, n.Variance); @@ -200,13 +200,13 @@ namespace MathNet.Numerics.UnitTests.DistributionTests [Row(10.0, 10.0, 0.31622776601683794197697302588502426416723164097476643)] [Row(10.0, 1.0, 3.1622776601683793319988935444327185337195551393252168)] [Row(10.0, Double.PositiveInfinity, 0.0)] - public void CanGetStdDev(double shape, double invScale, double sdev) + public void ValidateStdDev(double shape, double invScale, double sdev) { var n = new Gamma(shape, invScale); AssertHelpers.AlmostEqual(sdev, n.StdDev, 15); } - [Test] + [Test, Ignore("Depending on Special Functions")] [Row(0.0, 0.0, Double.PositiveInfinity)] [Row(1.0, 0.1, 3.3025850929940456285068402234265387271634735938763824)] [Row(1.0, 1.0, 1.0)] @@ -365,7 +365,7 @@ namespace MathNet.Numerics.UnitTests.DistributionTests var e = ied.Take(5).ToArray(); } - [Test] + [Test, Ignore("Depending on Special Functions")] [Row(0.0, 0.0, 0.0, 0.0)] [Row(0.0, 0.0, 1.0, 0.0)] [Row(0.0, 0.0, 10.0, 0.0)] diff --git a/src/Managed.UnitTests/Managed.UnitTests.csproj b/src/Managed.UnitTests/Managed.UnitTests.csproj index f140563f..589a731c 100644 --- a/src/Managed.UnitTests/Managed.UnitTests.csproj +++ b/src/Managed.UnitTests/Managed.UnitTests.csproj @@ -77,6 +77,7 @@ + diff --git a/src/Managed.UnitTests/SpecialFunctionsTest/SpecialFunctionsTests.cs b/src/Managed.UnitTests/SpecialFunctionsTest/SpecialFunctionsTests.cs new file mode 100644 index 00000000..2f5cf2e1 --- /dev/null +++ b/src/Managed.UnitTests/SpecialFunctionsTest/SpecialFunctionsTests.cs @@ -0,0 +1,81 @@ +// +// Math.NET Numerics, part of the Math.NET Project +// http://mathnet.opensourcedotnet.info +// +// Copyright (c) 2009 Math.NET +// +// Permission is hereby granted, free of charge, to any person +// obtaining a copy of this software and associated documentation +// files (the "Software"), to deal in the Software without +// restriction, including without limitation the rights to use, +// copy, modify, merge, publish, distribute, sublicense, and/or sell +// copies of the Software, and to permit persons to whom the +// Software is furnished to do so, subject to the following +// conditions: +// +// The above copyright notice and this permission notice shall be +// included in all copies or substantial portions of the Software. +// +// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, +// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES +// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND +// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT +// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, +// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING +// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR +// OTHER DEALINGS IN THE SOFTWARE. +// + +namespace MathNet.Numerics.UnitTests.SpecialFunctionTests +{ + using System; + using MbUnit.Framework; + using MathNet.Numerics; + + class SpecialFunctionsTests + { + [Test] + [Row(Double.NaN, Double.NaN)] + [Row(0.1, 2.2527126517342059020062379568954763844479865649307379)] + [Row(1.0, 0.0)] + [Row(1.5, -0.12078223763524522234551844578164721225185272790259947)] + [Row(Constants.Pi / 2, -0.11590380084550241329912089415904874214542604767006895)] + [Row(2.0, 0.0)] + [Row(2.5, 0.28468287047291915963249466968270192432013769555989498)] + [Row(3.0, 0.693147180559945309417232121458176568075500134360255)] + [Row(Constants.Pi, 0.82769459232343710152957855845235995115350173412073715)] + [Row(3.5, 1.2009736023470742248160218814507129957702389154681574)] + [Row(4.0, 1.7917594692280550008124773583807022727229906921830034)] + [Row(4.5, 2.4537365708424422205041425034357161573318235106897606)] + [Row(5.0, 3.1780538303479456196469416012970554088739909609035161)] + [Row(5.5, 3.9578139676187162938774008558225909985513044919750065)] + [Row(10.1, 13.02752673863323715481371189614224148681183971709386)] + public void GammaLn(double x, double f) + { + AssertHelpers.AlmostEqual(f, SpecialFunctions.GammaLn(x), 14); + } + + [Test] + [Row(Double.NaN, Double.NaN)] + [Row(-1.5, 2.3632718012073547030642233111215269103967326081631802)] + [Row(-0.5, -3.544907701811032054596334966682290365595098912244773)] + [Row(0.1, 9.5135076986687312858079798958252325009137161063903012)] + [Row(1.0, 1.0)] + [Row(1.5, 0.88622692545275801364908374167057259139877472806119326)] + [Row(Constants.Pi / 2, 0.89056089038153932801065963535912100593354196288475879)] + [Row(2.0, 1.0)] + [Row(2.5, 1.3293403881791370204736256125058588870981620920917912)] + [Row(3.0, 2.0)] + [Row(Constants.Pi, 2.2880377953400324179595889090602339228896881533562229)] + [Row(3.5, 3.3233509704478425511840640312646472177454052302294767)] + [Row(4.0, 6.0)] + [Row(4.5, 11.631728396567448929144224109426265262108918305803166)] + [Row(5.0, 24.0)] + [Row(5.5, 52.342777784553520181149008492418193679490132376114268)] + [Row(10.1, 454760.75144158558537612486797710217749925965322893332)] + public void Gamma(double x, double f) + { + AssertHelpers.AlmostEqual(f, SpecialFunctions.Gamma(x), 13); + } + } +} diff --git a/src/Managed/Constants.cs b/src/Managed/Constants.cs index df580dee..8677e005 100644 --- a/src/Managed/Constants.cs +++ b/src/Managed/Constants.cs @@ -98,6 +98,9 @@ namespace MathNet.Numerics /// The number log(sqrt(2*pi*e)) public const double LogSqrt2PiE = 1.4189385332046727417803297364056176398613974736378d; + /// The number log(2 * sqrt(e / pi)) + public const double LogTwoSqrtEOverPi = 0.6207822376352452223455184457816472122518527279025978; + /// The number 1/pi public const double InvPi = 0.31830988618379067153776752674502872406891929148091d; @@ -113,6 +116,9 @@ namespace MathNet.Numerics /// The number 2/sqrt(pi) public const double TwoInvSqrtPi = 1.1283791670955125738961589031215451716881012586580d; + /// The number 2 * sqrt(e / pi) + public const double TwoSqrtEOverPi = 1.8603827342052657173362492472666631120594218414085755; + /// The number (pi)/180 - factor to convert from Degree (deg) to Radians (rad). /// /// diff --git a/src/Managed/Distributions/Continuous/Beta.cs b/src/Managed/Distributions/Continuous/Beta.cs index c205bd8f..0d62e2a1 100644 --- a/src/Managed/Distributions/Continuous/Beta.cs +++ b/src/Managed/Distributions/Continuous/Beta.cs @@ -53,6 +53,11 @@ namespace MathNet.Numerics.Distributions /// private double _shapeB; + /// + /// The distribution's random number generator. + /// + private Random _random; + /// /// Initializes a new instance of the Beta distribution. /// @@ -81,7 +86,7 @@ namespace MathNet.Numerics.Distributions /// True when the parameters are valid, false otherwise. private static bool IsValidParameterSet(double a, double b) { - if (a < 0.0 || b < 0.0) + if (a < 0.0 || b < 0.0 || Double.IsNaN(a) || Double.IsNaN(b)) { return false; } @@ -129,14 +134,60 @@ namespace MathNet.Numerics.Distributions /// /// Gets or sets the random number generator which is used to draw random samples. /// - public Random RandomSource { get; set; } + public Random RandomSource + { + get + { + return _random; + } + + set + { + if (value == null) + { + throw new ArgumentNullException(); + } + + _random = value; + } + } /// /// Gets the mean of the Beta distribution. /// public double Mean { - get { return _shapeA / (_shapeA + _shapeB); } + get + { + if(_shapeA == 0.0 && _shapeB == 0.0) + { + return 0.5; + } + else if(_shapeA == 0.0) + { + return 0.0; + } + else if(_shapeB == 0.0) + { + return 1.0; + } + else if(Double.IsPositiveInfinity(_shapeA) && Double.IsPositiveInfinity(_shapeB)) + { + return 0.5; + } + else if (Double.IsPositiveInfinity(_shapeA)) + { + return 1.0; + } + else if (Double.IsPositiveInfinity(_shapeB)) + { + return 0.0; + } + else + { + return _shapeA / (_shapeA + _shapeB); + } + } } /// @@ -144,7 +195,10 @@ namespace MathNet.Numerics.Distributions /// public double Variance { - get { return (_shapeA * _shapeB) / ((_shapeA + _shapeB) * (_shapeA + _shapeB) * (_shapeA + _shapeB + 1.0)); } + get + { + return (_shapeA * _shapeB) / ((_shapeA + _shapeB) * (_shapeA + _shapeB) * (_shapeA + _shapeB + 1.0)); + } } /// @@ -185,11 +239,45 @@ namespace MathNet.Numerics.Distributions #region IContinuousDistribution implementation /// - /// Gets the mode of the Beta distribution. + /// Gets the mode of the Beta distribution; when there are multiple answers, this routine will return 0.5. /// public double Mode { - get { return (_shapeA - 1) / (_shapeA + _shapeB - 2); } + get + { + if (_shapeA == 0.0 && _shapeB == 0.0) + { + return 0.5; + } + else if (_shapeA == 0.0) + { + return 0.0; + } + else if (_shapeB == 0.0) + { + return 1.0; + } + else if (Double.IsPositiveInfinity(_shapeA) && Double.IsPositiveInfinity(_shapeB)) + { + return 0.5; + } + else if (Double.IsPositiveInfinity(_shapeA)) + { + return 1.0; + } + else if (Double.IsPositiveInfinity(_shapeB)) + { + return 0.0; + } + else if(_shapeA == 1.0 && _shapeB == 1.0) + { + return 0.5; + } + else + { + return (_shapeA - 1) / (_shapeA + _shapeB - 2); + } + } } /// diff --git a/src/Managed/Distributions/Continuous/ContinuousUniform.cs b/src/Managed/Distributions/Continuous/ContinuousUniform.cs index 6c167e54..a6584427 100644 --- a/src/Managed/Distributions/Continuous/ContinuousUniform.cs +++ b/src/Managed/Distributions/Continuous/ContinuousUniform.cs @@ -53,6 +53,11 @@ namespace MathNet.Numerics.Distributions /// private double _upper; + /// + /// The distribution's random number generator. + /// + private Random _random; + /// /// Initializes a new instance of the ContinuousUniform class with lower bound 0 and upper bound 1. /// @@ -156,7 +161,23 @@ namespace MathNet.Numerics.Distributions /// /// Gets or sets the random number generator which is used to draw random samples. /// - public Random RandomSource { get; set; } + public Random RandomSource + { + get + { + return _random; + } + + set + { + if (value == null) + { + throw new ArgumentNullException(); + } + + _random = value; + } + } /// /// Gets the mean of the distribution. diff --git a/src/Managed/Distributions/Continuous/Gamma.cs b/src/Managed/Distributions/Continuous/Gamma.cs index 379cb71b..1e66262c 100644 --- a/src/Managed/Distributions/Continuous/Gamma.cs +++ b/src/Managed/Distributions/Continuous/Gamma.cs @@ -60,6 +60,11 @@ namespace MathNet.Numerics.Distributions /// private double _invScale; + /// + /// The distribution's random number generator. + /// + private Random _random; + /// /// Initializes a new instance of the Gamma distribution. /// @@ -190,7 +195,23 @@ namespace MathNet.Numerics.Distributions /// /// Gets or sets the random number generator which is used to draw random samples. /// - public Random RandomSource { get; set; } + public Random RandomSource + { + get + { + return _random; + } + + set + { + if (value == null) + { + throw new ArgumentNullException(); + } + + _random = value; + } + } /// /// Gets the mean of the Gamma distribution. @@ -203,6 +224,10 @@ namespace MathNet.Numerics.Distributions { return _shape; } + else if(_invScale == 0.0 && _shape == 0.0) + { + return Double.NaN; + } else { return _shape / _invScale; diff --git a/src/Managed/Distributions/Continuous/Normal.cs b/src/Managed/Distributions/Continuous/Normal.cs index f7125140..1d2f15d6 100644 --- a/src/Managed/Distributions/Continuous/Normal.cs +++ b/src/Managed/Distributions/Continuous/Normal.cs @@ -53,6 +53,11 @@ namespace MathNet.Numerics.Distributions /// private double _stdDev; + /// + /// The distribution's random number generator. + /// + private Random _random; + /// /// Initializes a new instance of the Normal class. This is a normal distribution with mean 0.0 /// and standard deviation 1.0. The distribution will @@ -181,7 +186,23 @@ namespace MathNet.Numerics.Distributions /// /// Gets or sets the random number generator which is used to draw random samples. /// - public Random RandomSource { get; set; } + public Random RandomSource + { + get + { + return _random; + } + + set + { + if (value == null) + { + throw new ArgumentNullException(); + } + + _random = value; + } + } /// /// Gets or sets the mean of the normal distribution. diff --git a/src/Managed/SpecialFunctions.cs b/src/Managed/SpecialFunctions.cs index 76e1e74e..858d1947 100644 --- a/src/Managed/SpecialFunctions.cs +++ b/src/Managed/SpecialFunctions.cs @@ -36,6 +36,34 @@ namespace MathNet.Numerics /// public static partial class SpecialFunctions { + /// + /// The order of the GammaLn approximation. + /// + private const int Gamma_n = 10; + + /// + /// Auxiliary variable when evaluating the GammaLn function. + /// + private const double Gamma_r = 10.900511; + + /// + /// Polynomial coefficients for the GammaLn approximation. + /// + private static readonly double[] Gamma_dk = new double[] + { + 2.48574089138753565546e-5, + 1.05142378581721974210, + -3.45687097222016235469, + 4.51227709466894823700, + -2.98285225323576655721, + 1.05639711577126713077, + -1.95428773191645869583e-1, + 1.70970543404441224307e-2, + -5.71926117404305781283e-4, + 4.63399473359905636708e-6, + -2.71994908488607703910e-9 + }; + /// /// Computes the hypotenuse of a right angle triangle. /// @@ -59,36 +87,93 @@ namespace MathNet.Numerics return 0d; } - - public static double BetaLn(double a, double b) + /// + /// Computes the logarithm of the Gamma function. + /// + /// The argument of the gamma function. + /// The logarithm of the gamma function. + /// + /// This implementation of the computation of the gamma and logarithm of the gamma function follows the derivation in + /// "An Analysis Of The Lanczos Gamma Approximation", Glendon Ralph Pugh, 2004. + /// We use the implementation listed on p. 116 which achieves an accuracy of 16 floating point digits. Although 16 digit accuracy + /// should be sufficient for double values, improving accuracy is possible (see p. 126 in Pugh). + /// Our unit tests suggest that the accuracy of the Gamma function is correct up to 14 floating point digits. + /// + public static double GammaLn(double z) { - return Double.NaN; + if (z < 0.5) + { + double s = Gamma_dk[0]; + for (int i = 1; i <= Gamma_n; i++) + { + 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); + } + else + { + double s = Gamma_dk[0]; + for (int i = 1; i <= Gamma_n; i++) + { + 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); + } } - - public static double BetaRegularized(double a, double b, double x) + /// + /// Computes the Gamma function. + /// + /// The argument of the gamma function. + /// The logarithm of the gamma function. + /// + /// + /// This implementation of the computation of the gamma and logarithm of the gamma function follows the derivation in + /// "An Analysis Of The Lanczos Gamma Approximation", Glendon Ralph Pugh, 2004. + /// We use the implementation listed on p. 116 which should achieve an accuracy of 16 floating point digits. Although 16 digit accuracy + /// should be sufficient for double values, improving accuracy is possible (see p. 126 in Pugh). + /// + /// Our unit tests suggest that the accuracy of the Gamma function is correct up to 13 floating point digits. + /// + public static double Gamma(double z) { - return Double.NaN; + if (z < 0.5) + { + double s = Gamma_dk[0]; + for (int i = 1; i <= Gamma_n; i++) + { + 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 + { + double s = Gamma_dk[0]; + for (int i = 1; i <= Gamma_n; i++) + { + s += Gamma_dk[i] / (z + i - 1.0); + } + return s * Constants.TwoSqrtEOverPi * Math.Pow((z - 0.5 + Gamma_r) / Math.E, z - 0.5); + } } - public static double DiGamma(double x) + public static double IncompleteGamma(double x, double z, bool reg) { - return Double.NaN; + throw new NotImplementedException(); } - - public static double Gamma(double x) + public static double BetaLn(double a, double b) { - return Double.NaN; + throw new NotImplementedException(); } - public static double GammaLn(double x) + public static double BetaRegularized(double a, double b, double x) { - return Double.NaN; + throw new NotImplementedException(); } - public static double IncompleteGamma(double x, double z, bool reg) + public static double DiGamma(double x) { - return Double.NaN; + throw new NotImplementedException(); } } } diff --git a/src/Native.UnitTests/Native.UnitTests.csproj b/src/Native.UnitTests/Native.UnitTests.csproj index 214689f2..1d0c56ae 100644 --- a/src/Native.UnitTests/Native.UnitTests.csproj +++ b/src/Native.UnitTests/Native.UnitTests.csproj @@ -107,6 +107,9 @@ SpecialFunctionsTest\ErfTests.cs + + SpecialFunctionsTest\SpecialFunctionsTests.cs + StatisticsTests\DescriptiveStatisticsTests.cs