diff --git a/src/Numerics/Numerics.csproj b/src/Numerics/Numerics.csproj index c008f289..c46a500a 100644 --- a/src/Numerics/Numerics.csproj +++ b/src/Numerics/Numerics.csproj @@ -104,6 +104,7 @@ + @@ -412,7 +413,7 @@ - + diff --git a/src/Numerics/SpecialFunctions.cs b/src/Numerics/SpecialFunctions.cs deleted file mode 100644 index 5eeef137..00000000 --- a/src/Numerics/SpecialFunctions.cs +++ /dev/null @@ -1,205 +0,0 @@ -// -// Math.NET Numerics, part of the Math.NET Project -// http://numerics.mathdotnet.com -// http://github.com/mathnet/mathnet-numerics -// http://mathnetnumerics.codeplex.com -// -// Copyright (c) 2009-2011 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. -// - -// -// Cephes Math Library, Stephen L. Moshier -// ALGLIB 2.0.1, Sergey Bochkanov -// - -namespace MathNet.Numerics -{ - using System; - using Properties; - - /// - /// This class implements a collection of special function evaluations for double precision. This class - /// has a static constructor which will precompute a small number of values for faster runtime computations. - /// - public static partial class SpecialFunctions - { - /// - /// Initializes static members of the SpecialFunctions class. - /// - static SpecialFunctions() - { - InitializeFactorial(); - } - - /// - /// Computes the 'th Harmonic number. - /// - /// The Harmonic number which needs to be computed. - /// The t'th Harmonic number. - public static double Harmonic(int t) - { - return Constants.EulerMascheroni + DiGamma(t + 1.0); - } - - /// - /// Compute the generalized harmonic number of order n of m. (1 + 1/2^m + 1/3^m + ... + 1/n^m) - /// - /// The order parameter. - /// The power parameter. - /// General Harmonic number. - public static double GeneralHarmonic(int n, double m) - { - double sum = 0; - for (int i = 0; i < n; i++) - { - sum += Math.Pow(i + 1, -m); - } - return sum; - } - - /// - /// Computes the Digamma function which is mathematically defined as the derivative of the logarithm of the gamma function. - /// This implementation is based on - /// Jose Bernardo - /// Algorithm AS 103: - /// Psi ( Digamma ) Function, - /// Applied Statistics, - /// Volume 25, Number 3, 1976, pages 315-317. - /// Using the modifications as in Tom Minka's lightspeed toolbox. - /// - /// The argument of the digamma function. - /// The value of the DiGamma function at . - public static double DiGamma(double x) - { - const double C = 12.0; - const double D1 = -0.57721566490153286; - const double D2 = 1.6449340668482264365; - const double S = 1e-6; - const double S3 = 1.0 / 12.0; - const double S4 = 1.0 / 120.0; - const double S5 = 1.0 / 252.0; - const double S6 = 1.0 / 240.0; - const double S7 = 1.0 / 132.0; - - if (Double.IsNegativeInfinity(x) || Double.IsNaN(x)) - { - return Double.NaN; - } - - // Handle special cases. - if (x <= 0 && Math.Floor(x) == x) - { - return Double.NegativeInfinity; - } - - // Use inversion formula for negative numbers. - if (x < 0) - { - return DiGamma(1.0 - x) + (Math.PI / Math.Tan(-Math.PI * x)); - } - - if (x <= S) - { - return D1 - (1 / x) + (D2 * x); - } - - double result = 0; - while (x < C) - { - result -= 1 / x; - x++; - } - - if (x >= C) - { - var r = 1 / x; - result += Math.Log(x) - (0.5 * r); - r *= r; - - result -= r * (S3 - (r * (S4 - (r * (S5 - (r * (S6 - (r * S7)))))))); - } - - return result; - } - - /// - /// Computes the inverse Digamma function: this is the inverse of the logarithm of the gamma function. This function will - /// only return solutions that are positive. - /// This implementation is based on the bisection method. - /// - /// The argument of the inverse digamma function. - /// The positive solution to the inverse DiGamma function at . - public static double DiGammaInv(double p) - { - if (Double.IsNaN(p)) - { - return Double.NaN; - } - - if (Double.IsNegativeInfinity(p)) - { - return 0.0; - } - - if (Double.IsPositiveInfinity(p)) - { - return Double.PositiveInfinity; - } - - var x = Math.Exp(p); - for (var d = 1.0; d > 1.0e-15; d /= 2.0) - { - x += d * Math.Sign(p - DiGamma(x)); - } - - return x; - } - - /// - /// Computes the logit function. see: http://en.wikipedia.org/wiki/Logit - /// - /// The parameter for which to compute the logit function. This number should be - /// between 0 and 1. - /// The logarithm of divided by 1.0 - . - public static double Logit(double p) - { - if (p < 0.0 || p > 1.0) - { - throw new ArgumentOutOfRangeException(Resources.ArgumentBetween0And1); - } - - return Math.Log(p / (1.0 - p)); - } - - /// - /// Computes the logistic function. see: http://en.wikipedia.org/wiki/Logistic - /// - /// The parameter for which to compute the logistic function. - /// The logistic function of . - public static double Logistic(double p) - { - return 1.0 / (Math.Exp(-p) + 1.0); - } - } -} diff --git a/src/Numerics/SpecialFunctions/Beta.cs b/src/Numerics/SpecialFunctions/Beta.cs index cf7d668c..c96dfd93 100644 --- a/src/Numerics/SpecialFunctions/Beta.cs +++ b/src/Numerics/SpecialFunctions/Beta.cs @@ -33,7 +33,9 @@ // ALGLIB 2.0.1, Sergey Bochkanov // +// ReSharper disable CheckNamespace namespace MathNet.Numerics +// ReSharper restore CheckNamespace { using System; using Properties; diff --git a/src/Numerics/SpecialFunctions/Erf.cs b/src/Numerics/SpecialFunctions/Erf.cs index 4907f04b..4a62d00f 100644 --- a/src/Numerics/SpecialFunctions/Erf.cs +++ b/src/Numerics/SpecialFunctions/Erf.cs @@ -37,7 +37,9 @@ // * double ErfInvImpl(double p, double q, double s) // +// ReSharper disable CheckNamespace namespace MathNet.Numerics +// ReSharper restore CheckNamespace { using System; diff --git a/src/Numerics/SpecialFunctions/Factorial.cs b/src/Numerics/SpecialFunctions/Factorial.cs index 9308bc7f..fecdf339 100644 --- a/src/Numerics/SpecialFunctions/Factorial.cs +++ b/src/Numerics/SpecialFunctions/Factorial.cs @@ -28,7 +28,9 @@ // OTHER DEALINGS IN THE SOFTWARE. // +// ReSharper disable CheckNamespace namespace MathNet.Numerics +// ReSharper restore CheckNamespace { using System; using Properties; @@ -38,6 +40,14 @@ namespace MathNet.Numerics private const int FactorialMaxArgument = 170; private static double[] factorialCache; + /// + /// Initializes static members of the SpecialFunctions class. + /// + static SpecialFunctions() + { + InitializeFactorial(); + } + private static void InitializeFactorial() { factorialCache = new double[FactorialMaxArgument + 1]; diff --git a/src/Numerics/SpecialFunctions/Gamma.cs b/src/Numerics/SpecialFunctions/Gamma.cs index 97686982..39649d55 100644 --- a/src/Numerics/SpecialFunctions/Gamma.cs +++ b/src/Numerics/SpecialFunctions/Gamma.cs @@ -33,7 +33,9 @@ // ALGLIB 2.0.1, Sergey Bochkanov // +// ReSharper disable CheckNamespace namespace MathNet.Numerics +// ReSharper restore CheckNamespace { using System; @@ -362,5 +364,103 @@ namespace MathNet.Numerics return 1d - (Math.Exp(ax) * ans); } + + /// + /// Computes the Digamma function which is mathematically defined as the derivative of the logarithm of the gamma function. + /// This implementation is based on + /// Jose Bernardo + /// Algorithm AS 103: + /// Psi ( Digamma ) Function, + /// Applied Statistics, + /// Volume 25, Number 3, 1976, pages 315-317. + /// Using the modifications as in Tom Minka's lightspeed toolbox. + /// + /// The argument of the digamma function. + /// The value of the DiGamma function at . + public static double DiGamma(double x) + { + const double C = 12.0; + const double D1 = -0.57721566490153286; + const double D2 = 1.6449340668482264365; + const double S = 1e-6; + const double S3 = 1.0 / 12.0; + const double S4 = 1.0 / 120.0; + const double S5 = 1.0 / 252.0; + const double S6 = 1.0 / 240.0; + const double S7 = 1.0 / 132.0; + + if (Double.IsNegativeInfinity(x) || Double.IsNaN(x)) + { + return Double.NaN; + } + + // Handle special cases. + if (x <= 0 && Math.Floor(x) == x) + { + return Double.NegativeInfinity; + } + + // Use inversion formula for negative numbers. + if (x < 0) + { + return DiGamma(1.0 - x) + (Math.PI / Math.Tan(-Math.PI * x)); + } + + if (x <= S) + { + return D1 - (1 / x) + (D2 * x); + } + + double result = 0; + while (x < C) + { + result -= 1 / x; + x++; + } + + if (x >= C) + { + var r = 1 / x; + result += Math.Log(x) - (0.5 * r); + r *= r; + + result -= r * (S3 - (r * (S4 - (r * (S5 - (r * (S6 - (r * S7)))))))); + } + + return result; + } + + /// + /// Computes the inverse Digamma function: this is the inverse of the logarithm of the gamma function. This function will + /// only return solutions that are positive. + /// This implementation is based on the bisection method. + /// + /// The argument of the inverse digamma function. + /// The positive solution to the inverse DiGamma function at . + public static double DiGammaInv(double p) + { + if (Double.IsNaN(p)) + { + return Double.NaN; + } + + if (Double.IsNegativeInfinity(p)) + { + return 0.0; + } + + if (Double.IsPositiveInfinity(p)) + { + return Double.PositiveInfinity; + } + + var x = Math.Exp(p); + for (var d = 1.0; d > 1.0e-15; d /= 2.0) + { + x += d * Math.Sign(p - DiGamma(x)); + } + + return x; + } } } \ No newline at end of file diff --git a/src/Numerics/SpecialFunctions/Harmonic.cs b/src/Numerics/SpecialFunctions/Harmonic.cs new file mode 100644 index 00000000..2fe97214 --- /dev/null +++ b/src/Numerics/SpecialFunctions/Harmonic.cs @@ -0,0 +1,74 @@ +// +// Math.NET Numerics, part of the Math.NET Project +// http://numerics.mathdotnet.com +// http://github.com/mathnet/mathnet-numerics +// http://mathnetnumerics.codeplex.com +// +// Copyright (c) 2009-2011 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. +// + +// +// Cephes Math Library, Stephen L. Moshier +// ALGLIB 2.0.1, Sergey Bochkanov +// + +// ReSharper disable CheckNamespace +namespace MathNet.Numerics +// ReSharper restore CheckNamespace +{ + using System; + + /// + /// This partial implementation of the SpecialFunctions class contains all methods related to the harmonic function. + /// + public static partial class SpecialFunctions + { + + /// + /// Computes the 'th Harmonic number. + /// + /// The Harmonic number which needs to be computed. + /// The t'th Harmonic number. + public static double Harmonic(int t) + { + return Constants.EulerMascheroni + DiGamma(t + 1.0); + } + + /// + /// Compute the generalized harmonic number of order n of m. (1 + 1/2^m + 1/3^m + ... + 1/n^m) + /// + /// The order parameter. + /// The power parameter. + /// General Harmonic number. + public static double GeneralHarmonic(int n, double m) + { + double sum = 0; + for (int i = 0; i < n; i++) + { + sum += Math.Pow(i + 1, -m); + } + return sum; + } + } +} diff --git a/src/Numerics/SpecialFunctions/Logistic.cs b/src/Numerics/SpecialFunctions/Logistic.cs new file mode 100644 index 00000000..aae6ae6e --- /dev/null +++ b/src/Numerics/SpecialFunctions/Logistic.cs @@ -0,0 +1,74 @@ +// +// Math.NET Numerics, part of the Math.NET Project +// http://numerics.mathdotnet.com +// http://github.com/mathnet/mathnet-numerics +// http://mathnetnumerics.codeplex.com +// +// Copyright (c) 2009-2011 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. +// + +// +// Cephes Math Library, Stephen L. Moshier +// ALGLIB 2.0.1, Sergey Bochkanov +// + +// ReSharper disable CheckNamespace +namespace MathNet.Numerics +// ReSharper restore CheckNamespace +{ + using System; + using Properties; + + /// + /// This partial implementation of the SpecialFunctions class contains all methods related to the logistic function. + /// + public static partial class SpecialFunctions + { + /// + /// Computes the logistic function. see: http://en.wikipedia.org/wiki/Logistic + /// + /// The parameter for which to compute the logistic function. + /// The logistic function of . + public static double Logistic(double p) + { + return 1.0 / (Math.Exp(-p) + 1.0); + } + + /// + /// Computes the logit function, the inverse of the sigmoid logistic function. see: http://en.wikipedia.org/wiki/Logit + /// + /// The parameter for which to compute the logit function. This number should be + /// between 0 and 1. + /// The logarithm of divided by 1.0 - . + public static double Logit(double p) + { + if (p < 0.0 || p > 1.0) + { + throw new ArgumentOutOfRangeException(Resources.ArgumentBetween0And1); + } + + return Math.Log(p / (1.0 - p)); + } + } +} diff --git a/src/Numerics/SpecialFunctions/Stability.cs b/src/Numerics/SpecialFunctions/Stability.cs index 15a858ba..7cf72d85 100644 --- a/src/Numerics/SpecialFunctions/Stability.cs +++ b/src/Numerics/SpecialFunctions/Stability.cs @@ -28,7 +28,9 @@ // OTHER DEALINGS IN THE SOFTWARE. // +// ReSharper disable CheckNamespace namespace MathNet.Numerics +// ReSharper restore CheckNamespace { using System; using System.Numerics; diff --git a/src/Portable/Portable.csproj b/src/Portable/Portable.csproj index 98d20e2e..f44e3c78 100644 --- a/src/Portable/Portable.csproj +++ b/src/Portable/Portable.csproj @@ -984,9 +984,6 @@ Sorting.cs - - SpecialFunctions.cs - SpecialFunctions\Beta.cs @@ -999,6 +996,12 @@ SpecialFunctions\Gamma.cs + + SpecialFunctions\Harmonic.cs + + + SpecialFunctions\Logistic.cs + SpecialFunctions\Stability.cs