From f5f81a11ab472356a8c001361be0d9af81c5a595 Mon Sep 17 00:00:00 2001 From: Jurgen Van Gael Date: Thu, 27 Aug 2009 03:50:36 +0800 Subject: [PATCH] Fixed DiGammaInv bug for 0 and infinity arguments. Added more more special functions. Signed-off-by: jvangael Signed-off-by: jvangael --- src/Numerics/SpecialFunctions/Erf.cs | 43 ++++++++++++++++ src/Numerics/SpecialFunctions/Factorial.cs | 43 ++++++++++++++++ .../SpecialFunctionsTests/ErfTests.cs | 26 +++++++++- .../SpecialFunctionsTests/FactorialTest.cs | 9 ++++ .../SpecialFunctionsTests.cs | 49 ++++++++++++++++++- 5 files changed, 167 insertions(+), 3 deletions(-) diff --git a/src/Numerics/SpecialFunctions/Erf.cs b/src/Numerics/SpecialFunctions/Erf.cs index 3f626fad..5d48bea8 100644 --- a/src/Numerics/SpecialFunctions/Erf.cs +++ b/src/Numerics/SpecialFunctions/Erf.cs @@ -112,6 +112,49 @@ namespace MathNet.Numerics return ErfImp(x, true); } + ///Calculates the inverse error function evaluated at z. + /// The inverse error function evaluated at given value. + /// + /// + /// returns Double.PositiveInfinity if z >= 1.0. + /// returns Double.NegativeInfinity if z <= -1.0. + /// + /// + ///Calculates the inverse error function evaluated at z. + ///value to evaluate. + ///the inverse error function evaluated at Z. + public static double ErfInv(double z) + { + if (z == 0.0) + { + return 0.0; + } + + if (z >= 1.0) + { + return double.PositiveInfinity; + } + if (z <= -1.0) + { + return double.NegativeInfinity; + } + + double p, q, s; + if (z < 0) + { + p = -z; + q = 1 - p; + s = -1; + } + else + { + p = z; + q = 1 - z; + s = 1; + } + + return ErfInvImpl(p, q, s); + } /// /// Implementation of the error function. diff --git a/src/Numerics/SpecialFunctions/Factorial.cs b/src/Numerics/SpecialFunctions/Factorial.cs index 9d7c0e3d..e0f32638 100644 --- a/src/Numerics/SpecialFunctions/Factorial.cs +++ b/src/Numerics/SpecialFunctions/Factorial.cs @@ -29,6 +29,7 @@ namespace MathNet.Numerics { using System; + using Properties; public partial class SpecialFunctions { @@ -127,5 +128,47 @@ namespace MathNet.Numerics return FactorialLn(n) - FactorialLn(k) - FactorialLn(n - k); } + + /// + /// Computes the multinomial coefficient: n choose n1, n2, n3, ... + /// + /// A nonnegative value n. + /// An array of nonnegative values that sum to . + /// The multinomial coefficient. + /// if is . + /// If or any of the are negative. + /// If the sum of all is not equal to . + public static double Multinomial(int n, int[] ni) + { + if (n < 0) + { + throw new ArgumentException(Resources.ArgumentMustBePositive, "n"); + } + if (ni == null) + { + throw new ArgumentNullException("ni"); + } + + int sum = 0; + double ret = FactorialLn(n); + for (int i = 0; i < ni.Length; i++) + { + if (ni[i] < 0) + { + throw new ArgumentException(Resources.ArgumentMustBePositive, "ni[" + i + "]"); + } + + ret -= FactorialLn(ni[i]); + sum += ni[i]; + } + + // Before returning, check that the sum of all elements was equal to n. + if (sum != n) + { + throw new ArgumentException(Resources.ArgumentParameterSetInvalid , "ni"); + } + + return System.Math.Floor(0.5 + System.Math.Exp(ret)); + } } } \ No newline at end of file diff --git a/src/UnitTests/SpecialFunctionsTests/ErfTests.cs b/src/UnitTests/SpecialFunctionsTests/ErfTests.cs index d1ab0e05..9d0ed9f9 100644 --- a/src/UnitTests/SpecialFunctionsTests/ErfTests.cs +++ b/src/UnitTests/SpecialFunctionsTests/ErfTests.cs @@ -26,7 +26,7 @@ // OTHER DEALINGS IN THE SOFTWARE. // -namespace MathNet.Numerics.UnitTests.SpecialFunctionTests +namespace MathNet.Numerics.UnitTests.SpecialFunctionsTests { using MbUnit.Framework; using MathNet.Numerics; @@ -106,5 +106,29 @@ namespace MathNet.Numerics.UnitTests.SpecialFunctionTests { AssertHelpers.AlmostEqual(f, SpecialFunctions.ErfcInv(x), 8); } + + [Test] + [Row(double.NaN, double.NaN)] + [Row(-1.0, -0.84270079294971486934122063508260925929606699796630291)] + [Row(0.0, 0.0)] + [Row(1e-15, 0.0000000000000011283791670955126615773132947717431253912942469337536)] + [Row(0.1, 0.1124629160182848984047122510143040617233925185058162)] + [Row(0.2, 0.22270258921047846617645303120925671669511570710081967)] + [Row(0.3, 0.32862675945912741618961798531820303325847175931290341)] + [Row(0.4, 0.42839235504666847645410962730772853743532927705981257)] + [Row(0.5, 0.5204998778130465376827466538919645287364515757579637)] + [Row(1.0, 0.84270079294971486934122063508260925929606699796630291)] + [Row(1.5, 0.96610514647531072706697626164594785868141047925763678)] + [Row(2.0, 0.99532226501895273416206925636725292861089179704006008)] + [Row(2.5, 0.99959304798255504106043578426002508727965132259628658)] + [Row(3.0, 0.99997790950300141455862722387041767962015229291260075)] + [Row(4.0, 0.99999998458274209971998114784032651311595142785474641)] + [Row(5.0, 0.99999999999846254020557196514981165651461662110988195)] + [Row(double.PositiveInfinity, 1.0)] + [Row(double.NegativeInfinity, -1.0)] + public void ErfInv(double x, double f) + { + AssertHelpers.AlmostEqual(x, SpecialFunctions.ErfInv(f), 6); + } } } diff --git a/src/UnitTests/SpecialFunctionsTests/FactorialTest.cs b/src/UnitTests/SpecialFunctionsTests/FactorialTest.cs index 7c1ec6b9..2df887fd 100644 --- a/src/UnitTests/SpecialFunctionsTests/FactorialTest.cs +++ b/src/UnitTests/SpecialFunctionsTests/FactorialTest.cs @@ -102,5 +102,14 @@ namespace MathNet.Numerics.UnitTests.SpecialFunctionsTests AssertHelpers.AlmostEqual(Math.Log(0), SpecialFunctions.BinomialLn(5, 7), 14); AssertHelpers.AlmostEqual(Math.Log(0), SpecialFunctions.BinomialLn(5, -7), 14); } + + [Test] + public void CanComputeMultinomial() + { + AssertHelpers.AlmostEqual(1, SpecialFunctions.Multinomial(1, new int[] { 1, 0 }), 14); + AssertHelpers.AlmostEqual(10, SpecialFunctions.Multinomial(5, new int[] { 3, 2 }), 14); + AssertHelpers.AlmostEqual(10, SpecialFunctions.Multinomial(5, new int[] { 2, 3 }), 14); + AssertHelpers.AlmostEqual(35, SpecialFunctions.Multinomial(7, new int[] { 3, 4 }), 14); + } } } \ No newline at end of file diff --git a/src/UnitTests/SpecialFunctionsTests/SpecialFunctionsTests.cs b/src/UnitTests/SpecialFunctionsTests/SpecialFunctionsTests.cs index c5a0e186..3e0d1728 100644 --- a/src/UnitTests/SpecialFunctionsTests/SpecialFunctionsTests.cs +++ b/src/UnitTests/SpecialFunctionsTests/SpecialFunctionsTests.cs @@ -1,4 +1,4 @@ -// +// // Math.NET Numerics, part of the Math.NET Project // http://mathnet.opensourcedotnet.info // @@ -26,7 +26,7 @@ // OTHER DEALINGS IN THE SOFTWARE. // -namespace MathNet.Numerics.UnitTests.SpecialFunctionTests +namespace MathNet.Numerics.UnitTests.SpecialFunctionsTests { using System; using MbUnit.Framework; @@ -103,6 +103,7 @@ namespace MathNet.Numerics.UnitTests.SpecialFunctionTests [Test] [Row(Double.NaN, Double.NaN)] + [Row(0.0, Double.NegativeInfinity)] [Row(0.1, -10.423754940411076232100295314502760886768558023951363)] [Row(1.0, -0.57721566490153286060651209008240243104215933593992359)] [Row(1.5, 0.036489973978576520559023667001244432806840395339565888)] @@ -117,9 +118,53 @@ namespace MathNet.Numerics.UnitTests.SpecialFunctionTests [Row(5.0, 1.5061176684318004727268212432509309022911739973934097)] [Row(5.5, 1.6110931485817511237336268416044190359814435699427405)] [Row(10.1, 2.2622143570941481235561593642219403924532310597356171)] + [Row(Double.PositiveInfinity, Double.PositiveInfinity)] public void DiGammaInv(double x, double f) { AssertHelpers.AlmostEqual(x, SpecialFunctions.DiGammaInv(f), 13); } + + /// + /// Compute the t'th harmonic number using a loop. + /// + private double ExactHarmonic(int t) + { + double r = 0.0; + for (int i = 1; i <= t; i++) + { + r += 1.0 / i; + } + return r; + } + + [Test] + [Row(1)] + [Row(2)] + [Row(4)] + [Row(8)] + [Row(16)] + [Row(100)] + [Row(1000)] + [Row(10000)] + [Row(100000)] + [Row(1000000)] + public void Harmonic(int i) + { + AssertHelpers.AlmostEqual(ExactHarmonic(i), SpecialFunctions.Harmonic(i), 13); + } + + [Test] + public void BetaLn() + { + AssertHelpers.AlmostEqual(System.Math.Log(0.5), SpecialFunctions.BetaLn(1.0, 2.0), 14); + AssertHelpers.AlmostEqual(System.Math.Log(1.0), SpecialFunctions.BetaLn(1.0, 1.0), 14); + } + + [Test] + public void Beta() + { + AssertHelpers.AlmostEqual(0.5, SpecialFunctions.Beta(1.0, 2.0), 14); + AssertHelpers.AlmostEqual(1.0, SpecialFunctions.Beta(1.0, 1.0), 14); + } } } \ No newline at end of file