Browse Source

Added Gamma and GammaLn with unit tests.

Signed-off-by: Christoph Ruegg <git@cdrnet.ch>
pull/2/head
jvangael 17 years ago
committed by Christoph Ruegg
parent
commit
a8f0634c92
  1. 8
      src/Managed.UnitTests/AssertHelpers.cs
  2. 18
      src/Managed.UnitTests/DistributionTests/Continuous/BetaTests.cs
  3. 14
      src/Managed.UnitTests/DistributionTests/Continuous/GammaTests.cs
  4. 1
      src/Managed.UnitTests/Managed.UnitTests.csproj
  5. 81
      src/Managed.UnitTests/SpecialFunctionsTest/SpecialFunctionsTests.cs
  6. 6
      src/Managed/Constants.cs
  7. 100
      src/Managed/Distributions/Continuous/Beta.cs
  8. 23
      src/Managed/Distributions/Continuous/ContinuousUniform.cs
  9. 27
      src/Managed/Distributions/Continuous/Gamma.cs
  10. 23
      src/Managed/Distributions/Continuous/Normal.cs
  11. 115
      src/Managed/SpecialFunctions.cs
  12. 3
      src/Native.UnitTests/Native.UnitTests.csproj

8
src/Managed.UnitTests/AssertHelpers.cs

@ -38,13 +38,19 @@ namespace MathNet.Numerics.UnitTests
class AssertHelpers
{
/// <summary>
/// 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
/// <paramref name="expected"/> and <paramref name="actual"/> are NaN then no assert is thrown.
/// </summary>
/// <param name="expected">The expected value.</param>
/// <param name="actual">The actual value.</param>
/// <param name="decimalPlaces">The number of decimal places to agree on.</param>
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)
{

18
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<double>(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)]

14
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<double>(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<double>(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)]

1
src/Managed.UnitTests/Managed.UnitTests.csproj

@ -77,6 +77,7 @@
<Compile Include="Properties\AssemblyInfo.cs" />
<Compile Include="SortingTests.cs" />
<Compile Include="SpecialFunctionsTest\ErfTests.cs" />
<Compile Include="SpecialFunctionsTest\SpecialFunctionsTests.cs" />
<Compile Include="StatisticsTests\DescriptiveStatisticsTests.cs" />
<Compile Include="StatisticsTests\StatisticsTests.cs" />
<Compile Include="StatisticsTests\StatTestData.cs" />

81
src/Managed.UnitTests/SpecialFunctionsTest/SpecialFunctionsTests.cs

@ -0,0 +1,81 @@
// <copyright file="SpecialFunctionsTests.cs" company="Math.NET">
// 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.
// </copyright>
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);
}
}
}

6
src/Managed/Constants.cs

@ -98,6 +98,9 @@ namespace MathNet.Numerics
/// <summary>The number log(sqrt(2*pi*e))</summary>
public const double LogSqrt2PiE = 1.4189385332046727417803297364056176398613974736378d;
/// <summary>The number log(2 * sqrt(e / pi))</summary>
public const double LogTwoSqrtEOverPi = 0.6207822376352452223455184457816472122518527279025978;
/// <summary>The number 1/pi</summary>
public const double InvPi = 0.31830988618379067153776752674502872406891929148091d;
@ -113,6 +116,9 @@ namespace MathNet.Numerics
/// <summary>The number 2/sqrt(pi)</summary>
public const double TwoInvSqrtPi = 1.1283791670955125738961589031215451716881012586580d;
/// <summary>The number 2 * sqrt(e / pi)</summary>
public const double TwoSqrtEOverPi = 1.8603827342052657173362492472666631120594218414085755;
/// <summary>The number (pi)/180 - factor to convert from Degree (deg) to Radians (rad).</summary>
/// <seealso cref="Trig.DegreeToRadian"/>
/// <seealso cref="Trig.RadianToDegree"/>

100
src/Managed/Distributions/Continuous/Beta.cs

@ -53,6 +53,11 @@ namespace MathNet.Numerics.Distributions
/// </summary>
private double _shapeB;
/// <summary>
/// The distribution's random number generator.
/// </summary>
private Random _random;
/// <summary>
/// Initializes a new instance of the Beta distribution.
/// </summary>
@ -81,7 +86,7 @@ namespace MathNet.Numerics.Distributions
/// <returns>True when the parameters are valid, false otherwise.</returns>
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
/// <summary>
/// Gets or sets the random number generator which is used to draw random samples.
/// </summary>
public Random RandomSource { get; set; }
public Random RandomSource
{
get
{
return _random;
}
set
{
if (value == null)
{
throw new ArgumentNullException();
}
_random = value;
}
}
/// <summary>
/// Gets the mean of the Beta distribution.
/// </summary>
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);
}
}
}
/// <summary>
@ -144,7 +195,10 @@ namespace MathNet.Numerics.Distributions
/// </summary>
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));
}
}
/// <summary>
@ -185,11 +239,45 @@ namespace MathNet.Numerics.Distributions
#region IContinuousDistribution implementation
/// <summary>
/// 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.
/// </summary>
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);
}
}
}
/// <summary>

23
src/Managed/Distributions/Continuous/ContinuousUniform.cs

@ -53,6 +53,11 @@ namespace MathNet.Numerics.Distributions
/// </summary>
private double _upper;
/// <summary>
/// The distribution's random number generator.
/// </summary>
private Random _random;
/// <summary>
/// Initializes a new instance of the ContinuousUniform class with lower bound 0 and upper bound 1.
/// </summary>
@ -156,7 +161,23 @@ namespace MathNet.Numerics.Distributions
/// <summary>
/// Gets or sets the random number generator which is used to draw random samples.
/// </summary>
public Random RandomSource { get; set; }
public Random RandomSource
{
get
{
return _random;
}
set
{
if (value == null)
{
throw new ArgumentNullException();
}
_random = value;
}
}
/// <summary>
/// Gets the mean of the distribution.

27
src/Managed/Distributions/Continuous/Gamma.cs

@ -60,6 +60,11 @@ namespace MathNet.Numerics.Distributions
/// </summary>
private double _invScale;
/// <summary>
/// The distribution's random number generator.
/// </summary>
private Random _random;
/// <summary>
/// Initializes a new instance of the Gamma distribution.
/// </summary>
@ -190,7 +195,23 @@ namespace MathNet.Numerics.Distributions
/// <summary>
/// Gets or sets the random number generator which is used to draw random samples.
/// </summary>
public Random RandomSource { get; set; }
public Random RandomSource
{
get
{
return _random;
}
set
{
if (value == null)
{
throw new ArgumentNullException();
}
_random = value;
}
}
/// <summary>
/// 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;

23
src/Managed/Distributions/Continuous/Normal.cs

@ -53,6 +53,11 @@ namespace MathNet.Numerics.Distributions
/// </summary>
private double _stdDev;
/// <summary>
/// The distribution's random number generator.
/// </summary>
private Random _random;
/// <summary>
/// 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
/// <summary>
/// Gets or sets the random number generator which is used to draw random samples.
/// </summary>
public Random RandomSource { get; set; }
public Random RandomSource
{
get
{
return _random;
}
set
{
if (value == null)
{
throw new ArgumentNullException();
}
_random = value;
}
}
/// <summary>
/// Gets or sets the mean of the normal distribution.

115
src/Managed/SpecialFunctions.cs

@ -36,6 +36,34 @@ namespace MathNet.Numerics
/// </summary>
public static partial class SpecialFunctions
{
/// <summary>
/// The order of the GammaLn approximation.
/// </summary>
private const int Gamma_n = 10;
/// <summary>
/// Auxiliary variable when evaluating the GammaLn function.
/// </summary>
private const double Gamma_r = 10.900511;
/// <summary>
/// Polynomial coefficients for the GammaLn approximation.
/// </summary>
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
};
/// <summary>
/// Computes the hypotenuse of a right angle triangle.
/// </summary>
@ -59,36 +87,93 @@ namespace MathNet.Numerics
return 0d;
}
public static double BetaLn(double a, double b)
/// <summary>
/// Computes the logarithm of the Gamma function.
/// </summary>
/// <param name="z">The argument of the gamma function.</param>
/// <returns>The logarithm of the gamma function.</returns>
/// <remarks>
/// <para>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).</para>
/// <para>Our unit tests suggest that the accuracy of the Gamma function is correct up to 14 floating point digits.</para>
/// </remarks>
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)
/// <summary>
/// Computes the Gamma function.
/// </summary>
/// <param name="z">The argument of the gamma function.</param>
/// <returns>The logarithm of the gamma function.</returns>
/// <remarks>
/// <para>
/// 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).
/// </para>
/// <para>Our unit tests suggest that the accuracy of the Gamma function is correct up to 13 floating point digits.</para>
/// </remarks>
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();
}
}
}

3
src/Native.UnitTests/Native.UnitTests.csproj

@ -107,6 +107,9 @@
<Compile Include="..\Managed.UnitTests\SpecialFunctionsTest\ErfTests.cs">
<Link>SpecialFunctionsTest\ErfTests.cs</Link>
</Compile>
<Compile Include="..\Managed.UnitTests\SpecialFunctionsTest\SpecialFunctionsTests.cs">
<Link>SpecialFunctionsTest\SpecialFunctionsTests.cs</Link>
</Compile>
<Compile Include="..\Managed.UnitTests\StatisticsTests\DescriptiveStatisticsTests.cs">
<Link>StatisticsTests\DescriptiveStatisticsTests.cs</Link>
</Compile>

Loading…
Cancel
Save