diff --git a/src/Numerics/Numerics.csproj b/src/Numerics/Numerics.csproj index e1153e2b..0823b1f9 100644 --- a/src/Numerics/Numerics.csproj +++ b/src/Numerics/Numerics.csproj @@ -456,6 +456,7 @@ + diff --git a/src/Numerics/SpecialFunctions.cs b/src/Numerics/SpecialFunctions.cs index d70dfab2..7fa3239a 100644 --- a/src/Numerics/SpecialFunctions.cs +++ b/src/Numerics/SpecialFunctions.cs @@ -3,7 +3,9 @@ // http://numerics.mathdotnet.com // http://github.com/mathnet/mathnet-numerics // http://mathnetnumerics.codeplex.com -// Copyright (c) 2009-2010 Math.NET +// +// 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 @@ -12,8 +14,10 @@ // 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 @@ -74,40 +78,6 @@ namespace MathNet.Numerics return sum; } - /// - /// Computes the logarithm of the Euler Beta function. - /// - /// The first Beta parameter, a positive real number. - /// The second Beta parameter, a positive real number. - /// The logarithm of the Euler Beta function evaluated at z,w. - /// If or are not positive. - public static double BetaLn(double z, double w) - { - if (z <= 0.0) - { - throw new ArgumentException(Resources.ArgumentMustBePositive, "z"); - } - - if (w <= 0.0) - { - throw new ArgumentException(Resources.ArgumentMustBePositive, "w"); - } - - return GammaLn(z) + GammaLn(w) - GammaLn(z + w); - } - - /// - /// Computes the Euler Beta function. - /// - /// The first Beta parameter, a positive real number. - /// The second Beta parameter, a positive real number. - /// The Euler Beta function evaluated at z,w. - /// If or are not positive. - public static double Beta(double z, double w) - { - return Math.Exp(BetaLn(z, w)); - } - /// /// Computes the Digamma function which is mathematically defined as the derivative of the logarithm of the gamma function. /// This implementation is based on @@ -206,128 +176,6 @@ namespace MathNet.Numerics return x; } - /// - /// Returns the lower incomplete (unregularized) beta function - /// I_x(a,b) = int(t^(a-1)*(1-t)^(b-1),t=0..x) for real a > 0, b > 0, 1 >= x >= 0. - /// - /// The first Beta parameter, a positive real number. - /// The second Beta parameter, a positive real number. - /// The upper limit of the integral. - /// The lower incomplete (unregularized) beta function. - public static double BetaIncomplete(double a, double b, double x) - { - return BetaRegularized(a, b, x) * Beta(a, b); - } - - /// - /// Returns the regularized lower incomplete beta function - /// I_x(a,b) = 1/Beta(a,b) * int(t^(a-1)*(1-t)^(b-1),t=0..x) for real a > 0, b > 0, 1 >= x >= 0. - /// - /// The first Beta parameter, a positive real number. - /// The second Beta parameter, a positive real number. - /// The upper limit of the integral. - /// The regularized lower incomplete beta function. - public static double BetaRegularized(double a, double b, double x) - { - if (a < 0.0) - { - throw new ArgumentOutOfRangeException("a", Resources.ArgumentNotNegative); - } - - if (b < 0.0) - { - throw new ArgumentOutOfRangeException("b", Resources.ArgumentNotNegative); - } - - if (x < 0.0 || x > 1.0) - { - throw new ArgumentOutOfRangeException("x", Resources.ArgumentInIntervalXYInclusive); - } - - var bt = (x == 0.0 || x == 1.0) - ? 0.0 - : Math.Exp(GammaLn(a + b) - GammaLn(a) - GammaLn(b) + (a * Math.Log(x)) + (b * Math.Log(1.0 - x))); - - var symmetryTransformation = x >= (a + 1.0) / (a + b + 2.0); - - /* Continued fraction representation */ - const int MaxIterations = 100; - var eps = Precision.DoubleMachinePrecision; - var fpmin = 0.0.Increment() / eps; - - if (symmetryTransformation) - { - x = 1.0 - x; - var swap = a; - a = b; - b = swap; - } - - var qab = a + b; - var qap = a + 1.0; - var qam = a - 1.0; - var c = 1.0; - var d = 1.0 - (qab * x / qap); - - if (Math.Abs(d) < fpmin) - { - d = fpmin; - } - - d = 1.0 / d; - var h = d; - - for (int m = 1, m2 = 2; m <= MaxIterations; m++, m2 += 2) - { - var aa = m * (b - m) * x / ((qam + m2) * (a + m2)); - d = 1.0 + (aa * d); - - if (Math.Abs(d) < fpmin) - { - d = fpmin; - } - - c = 1.0 + (aa / c); - if (Math.Abs(c) < fpmin) - { - c = fpmin; - } - - d = 1.0 / d; - h *= d * c; - aa = -(a + m) * (qab + m) * x / ((a + m2) * (qap + m2)); - d = 1.0 + (aa * d); - - if (Math.Abs(d) < fpmin) - { - d = fpmin; - } - - c = 1.0 + (aa / c); - - if (Math.Abs(c) < fpmin) - { - c = fpmin; - } - - d = 1.0 / d; - var del = d * c; - h *= del; - - if (Math.Abs(del - 1.0) <= eps) - { - if (symmetryTransformation) - { - return 1.0 - (bt * h / a); - } - - return bt * h / a; - } - } - - throw new ArgumentException(Resources.ArgumentTooLargeForIterationLimit); - } - /// /// Computes the logit function. see: http://en.wikipedia.org/wiki/Logit /// diff --git a/src/Numerics/SpecialFunctions/Beta.cs b/src/Numerics/SpecialFunctions/Beta.cs new file mode 100644 index 00000000..8b2dbe15 --- /dev/null +++ b/src/Numerics/SpecialFunctions/Beta.cs @@ -0,0 +1,199 @@ +// +// 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, Sergey Bochkanov +// + +namespace MathNet.Numerics +{ + using System; + using Properties; + + public static partial class SpecialFunctions + { + /// + /// Computes the logarithm of the Euler Beta function. + /// + /// The first Beta parameter, a positive real number. + /// The second Beta parameter, a positive real number. + /// The logarithm of the Euler Beta function evaluated at z,w. + /// If or are not positive. + public static double BetaLn(double z, double w) + { + if (z <= 0.0) + { + throw new ArgumentException(Resources.ArgumentMustBePositive, "z"); + } + + if (w <= 0.0) + { + throw new ArgumentException(Resources.ArgumentMustBePositive, "w"); + } + + return GammaLn(z) + GammaLn(w) - GammaLn(z + w); + } + + /// + /// Computes the Euler Beta function. + /// + /// The first Beta parameter, a positive real number. + /// The second Beta parameter, a positive real number. + /// The Euler Beta function evaluated at z,w. + /// If or are not positive. + public static double Beta(double z, double w) + { + return Math.Exp(BetaLn(z, w)); + } + + /// + /// Returns the lower incomplete (unregularized) beta function + /// I_x(a,b) = int(t^(a-1)*(1-t)^(b-1),t=0..x) for real a > 0, b > 0, 1 >= x >= 0. + /// + /// The first Beta parameter, a positive real number. + /// The second Beta parameter, a positive real number. + /// The upper limit of the integral. + /// The lower incomplete (unregularized) beta function. + public static double BetaIncomplete(double a, double b, double x) + { + return BetaRegularized(a, b, x) * Beta(a, b); + } + + /// + /// Returns the regularized lower incomplete beta function + /// I_x(a,b) = 1/Beta(a,b) * int(t^(a-1)*(1-t)^(b-1),t=0..x) for real a > 0, b > 0, 1 >= x >= 0. + /// + /// The first Beta parameter, a positive real number. + /// The second Beta parameter, a positive real number. + /// The upper limit of the integral. + /// The regularized lower incomplete beta function. + public static double BetaRegularized(double a, double b, double x) + { + if (a < 0.0) + { + throw new ArgumentOutOfRangeException("a", Resources.ArgumentNotNegative); + } + + if (b < 0.0) + { + throw new ArgumentOutOfRangeException("b", Resources.ArgumentNotNegative); + } + + if (x < 0.0 || x > 1.0) + { + throw new ArgumentOutOfRangeException("x", Resources.ArgumentInIntervalXYInclusive); + } + + var bt = (x == 0.0 || x == 1.0) + ? 0.0 + : Math.Exp(GammaLn(a + b) - GammaLn(a) - GammaLn(b) + (a * Math.Log(x)) + (b * Math.Log(1.0 - x))); + + var symmetryTransformation = x >= (a + 1.0) / (a + b + 2.0); + + /* Continued fraction representation */ + const int MaxIterations = 100; + var eps = Precision.DoubleMachinePrecision; + var fpmin = 0.0.Increment() / eps; + + if (symmetryTransformation) + { + x = 1.0 - x; + var swap = a; + a = b; + b = swap; + } + + var qab = a + b; + var qap = a + 1.0; + var qam = a - 1.0; + var c = 1.0; + var d = 1.0 - (qab * x / qap); + + if (Math.Abs(d) < fpmin) + { + d = fpmin; + } + + d = 1.0 / d; + var h = d; + + for (int m = 1, m2 = 2; m <= MaxIterations; m++, m2 += 2) + { + var aa = m * (b - m) * x / ((qam + m2) * (a + m2)); + d = 1.0 + (aa * d); + + if (Math.Abs(d) < fpmin) + { + d = fpmin; + } + + c = 1.0 + (aa / c); + if (Math.Abs(c) < fpmin) + { + c = fpmin; + } + + d = 1.0 / d; + h *= d * c; + aa = -(a + m) * (qab + m) * x / ((a + m2) * (qap + m2)); + d = 1.0 + (aa * d); + + if (Math.Abs(d) < fpmin) + { + d = fpmin; + } + + c = 1.0 + (aa / c); + + if (Math.Abs(c) < fpmin) + { + c = fpmin; + } + + d = 1.0 / d; + var del = d * c; + h *= del; + + if (Math.Abs(del - 1.0) <= eps) + { + if (symmetryTransformation) + { + return 1.0 - (bt * h / a); + } + + return bt * h / a; + } + } + + throw new ArgumentException(Resources.ArgumentTooLargeForIterationLimit); + } + } +} diff --git a/src/Numerics/SpecialFunctions/Gamma.cs b/src/Numerics/SpecialFunctions/Gamma.cs index 715711da..10a9f11a 100644 --- a/src/Numerics/SpecialFunctions/Gamma.cs +++ b/src/Numerics/SpecialFunctions/Gamma.cs @@ -36,7 +36,6 @@ namespace MathNet.Numerics { using System; - using Properties; public static partial class SpecialFunctions { diff --git a/src/Silverlight/Silverlight.csproj b/src/Silverlight/Silverlight.csproj index bc6446c8..cc8ba5f5 100644 --- a/src/Silverlight/Silverlight.csproj +++ b/src/Silverlight/Silverlight.csproj @@ -948,6 +948,9 @@ SpecialFunctions.cs + + SpecialFunctions\Beta.cs + SpecialFunctions\Erf.cs