diff --git a/src/Numerics/Distributions/Continuous/Weibull.cs b/src/Numerics/Distributions/Continuous/Weibull.cs index 50d8eb67..9618513e 100644 --- a/src/Numerics/Distributions/Continuous/Weibull.cs +++ b/src/Numerics/Distributions/Continuous/Weibull.cs @@ -55,6 +55,15 @@ namespace MathNet.Numerics.Distributions /// private double _scale; + /// + /// Reusable intermediate result 1 / ( ^ ) + /// + /// + /// By caching this parameter we can get slightly better numerics precision + /// in certain constellations without any additional computations. + /// + private double _scalePowShapeInv; + /// /// The distribution's random number generator. /// @@ -111,6 +120,7 @@ namespace MathNet.Numerics.Distributions _shape = shape; _scale = scale; + _scalePowShapeInv = Math.Pow(scale, -shape); } /// @@ -239,14 +249,12 @@ namespace MathNet.Numerics.Distributions { get { - if (_shape > 1.0) - { - return _scale * Math.Pow((_shape - 1.0) / _shape, 1.0 / _shape); - } - else + if (_shape <= 1.0) { return 0.0; } + + return _scale * Math.Pow((_shape - 1.0) / _shape, 1.0 / _shape); } } @@ -290,10 +298,8 @@ namespace MathNet.Numerics.Distributions { return _shape / _scale; } - else - { - return _shape * Math.Pow(x / _scale, _shape - 1.0) * Math.Exp(-Math.Pow(x / _scale, _shape)) / _scale; - } + + return _shape * Math.Pow(x / _scale, _shape - 1.0) * Math.Exp(-Math.Pow(x, _shape) * _scalePowShapeInv) / _scale; } return 0.0; @@ -312,10 +318,8 @@ namespace MathNet.Numerics.Distributions { return Math.Log(_shape) - Math.Log(_scale); } - else - { - return Math.Log(_shape) + (_shape - 1.0) * Math.Log(x / _scale) - Math.Pow(x / _scale, _shape) - Math.Log(_scale); - } + + return Math.Log(_shape) + (_shape - 1.0) * Math.Log(x / _scale) - (Math.Pow(x, _shape) * _scalePowShapeInv) - Math.Log(_scale); } return double.NegativeInfinity; @@ -328,12 +332,12 @@ namespace MathNet.Numerics.Distributions /// the cumulative density at . public double CumulativeDistribution(double x) { - if (x >= 0.0) + if (x < 0.0) { - return 1.0 - Math.Exp(-Math.Pow(x / _scale, _shape)); + return 0.0; } - return 0.0; + return -SpecialFunctions.ExponentialMinusOne(-Math.Pow(x, _shape) * _scalePowShapeInv); } /// diff --git a/src/Numerics/Numerics.csproj b/src/Numerics/Numerics.csproj index fee8bf91..26448f4b 100644 --- a/src/Numerics/Numerics.csproj +++ b/src/Numerics/Numerics.csproj @@ -112,6 +112,7 @@ + diff --git a/src/Numerics/SpecialFunctions.cs b/src/Numerics/SpecialFunctions.cs index 2b5a847c..1bb90334 100644 --- a/src/Numerics/SpecialFunctions.cs +++ b/src/Numerics/SpecialFunctions.cs @@ -73,29 +73,6 @@ namespace MathNet.Numerics InitializeFactorial(); } - /// - /// Computes the hypotenuse of a right angle triangle. - /// - /// The length of side a of the triangle. - /// The length of side b of the triangle. - /// Returns sqrt(a2 + b2) without underflow/overflow. - public static double Hypotenuse(double a, double b) - { - if (Math.Abs(a) > Math.Abs(b)) - { - double r = b / a; - return Math.Abs(a) * Math.Sqrt(1 + (r * r)); - } - - if (!b.AlmostZero()) - { - double r = a / b; - return Math.Abs(b) * Math.Sqrt(1 + (r * r)); - } - - return 0d; - } - /// /// Computes the logarithm of the Gamma function. /// diff --git a/src/Numerics/SpecialFunctions/Stability.cs b/src/Numerics/SpecialFunctions/Stability.cs new file mode 100644 index 00000000..7bd0acaa --- /dev/null +++ b/src/Numerics/SpecialFunctions/Stability.cs @@ -0,0 +1,120 @@ +// +// 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 +{ + using System; + + public partial class SpecialFunctions + { + /// + /// Numerically stable exponential minus one, i.e. x -> exp(x)-1 + /// + /// A number specifying a power. + /// Returns exp(power)-1. + public static double ExponentialMinusOne(double power) + { + double x = Math.Abs(power); + if (x > 0.1) + { + return Math.Exp(power) - 1.0; + } + + if (x < Precision.DoubleMachinePrecision) + { + return x; + } + + // Series Expansion to x^k / k! + int k = 0; + double term = 1.0; + return Series( + () => + { + k++; + term *= power; + term /= k; + return term; + } + ); + } + + /// + /// Numerically stable hypotenuse of a right angle triangle, i.e. (a,b) -> sqrt(a^2 + b^2) + /// + /// The length of side a of the triangle. + /// The length of side b of the triangle. + /// Returns sqrt(a2 + b2) without underflow/overflow. + public static double Hypotenuse(double a, double b) + { + if (Math.Abs(a) > Math.Abs(b)) + { + double r = b / a; + return Math.Abs(a) * Math.Sqrt(1 + (r * r)); + } + + if (b != 0.0) + { + // NOTE (ruegg): not "!b.AlmostZero()" to avoid convergence issues (e.g. in SVD algorithm) + double r = a / b; + return Math.Abs(b) * Math.Sqrt(1 + (r * r)); + } + + return 0d; + } + + /// + /// Numerically stable series summation + /// + /// provides the summands sequentially + /// Sum + private static double Series(Func nextSummand) + { + double compensation = 0.0; + double current; + double factor = 1 << 16; + + double sum = nextSummand(); + + do + { + // Kahan Summation + // NOTE (ruegg): do NOT optimize. Now, how to tell that the compiler? + current = nextSummand(); + double y = current - compensation; + double t = sum + y; + compensation = t - sum; + compensation -= y; + sum = t; + } + while (Math.Abs(sum) < Math.Abs(factor * current)); + + return sum; + } + } +} \ No newline at end of file diff --git a/src/UnitTests/DistributionTests/Continuous/WeibullTests.cs b/src/UnitTests/DistributionTests/Continuous/WeibullTests.cs index 5c56adef..23c5b578 100644 --- a/src/UnitTests/DistributionTests/Continuous/WeibullTests.cs +++ b/src/UnitTests/DistributionTests/Continuous/WeibullTests.cs @@ -282,7 +282,7 @@ namespace MathNet.Numerics.UnitTests.DistributionTests var e = ied.Take(5).ToArray(); } - [Test, Ignore("Catastrophic cancellation in one case. Fix this.")] + [Test] [Row(1.0, 0.1, 0.0, 0.0)] [Row(1.0, 0.1, 1.0, 0.99995460007023751514846440848443944938976208191113)] [Row(1.0, 0.1, 10.0, 0.99999999999999999999999999999999999999999996279924)]