functions: stable x -> exp(x)-1, used in weibull cdf

Signed-off-by: Christoph Ruegg <git@cdrnet.ch>

src/Numerics/Distributions/Continuous/Weibull.cs

@@ -55,6 +55,15 @@ namespace MathNet.Numerics.Distributions
         /// </summary>
         private double _scale;
 
+        /// <summary>
+        /// Reusable intermediate result 1 / (<see cref="_scale"/> ^ <see cref="_shape"/>)
+        /// </summary>
+        /// <remarks>
+        /// By caching this parameter we can get slightly better numerics precision
+        /// in certain constellations without any additional computations.
+        /// </remarks>
+        private double _scalePowShapeInv;
+
         /// <summary>
         /// The distribution's random number generator.
         /// </summary>
@@ -111,6 +120,7 @@ namespace MathNet.Numerics.Distributions
 
             _shape = shape;
             _scale = scale;
+            _scalePowShapeInv = Math.Pow(scale, -shape);
         }
 
         /// <summary>
@@ -239,14 +249,12 @@ namespace MathNet.Numerics.Distributions
         {
             get
             {
-                if (_shape > 1.0)
+                if (_shape <= 1.0)
-                {
-                    return _scale * Math.Pow((_shape - 1.0) / _shape, 1.0 / _shape);
-                }
-                else
                 {
                     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, _shape) * _scalePowShapeInv) / _scale;
-                    return _shape * Math.Pow(x / _scale, _shape - 1.0) * Math.Exp(-Math.Pow(x / _scale, _shape)) / _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, _shape) * _scalePowShapeInv) - Math.Log(_scale);
-                    return Math.Log(_shape) + (_shape - 1.0) * Math.Log(x / _scale) - Math.Pow(x / _scale, _shape) - Math.Log(_scale);
-                }
             }
 
             return double.NegativeInfinity;
@@ -328,12 +332,12 @@ namespace MathNet.Numerics.Distributions
         /// <returns>the cumulative density at <paramref name="x"/>.</returns>
         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);
         }
 
         /// <summary>

src/Numerics/Numerics.csproj

@@ -112,6 +112,7 @@
     <Compile Include="SpecialFunctions.cs" />
     <Compile Include="SpecialFunctions\Erf.cs" />
     <Compile Include="SpecialFunctions\Factorial.cs" />
+    <Compile Include="SpecialFunctions\Stability.cs" />
     <Compile Include="Statistics\DescriptiveStatistics.cs" />
     <Compile Include="Statistics\Statistics.cs" />
     <Compile Include="Threading\AggregateException.cs" />

src/Numerics/SpecialFunctions.cs

@@ -73,29 +73,6 @@ namespace MathNet.Numerics
             InitializeFactorial();
         }
 
-        /// <summary>
-        /// Computes the hypotenuse of a right angle triangle.
-        /// </summary>
-        /// <param name="a">The length of side a of the triangle.</param>
-        /// <param name="b">The length of side b of the triangle.</param>
-        /// <returns>Returns <code>sqrt(a<sup>2</sup> + b<sup>2</sup>)</code> without underflow/overflow.</returns>
-        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;
-        }
-
         /// <summary>
         /// Computes the logarithm of the Gamma function. 
         /// </summary>

src/Numerics/SpecialFunctions/Stability.cs

@@ -0,0 +1,120 @@
+// <copyright file="Stability.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
+{
+    using System;
+
+    public partial class SpecialFunctions
+    {
+        /// <summary>
+        /// Numerically stable exponential minus one, i.e. <code>x -> exp(x)-1</code>
+        /// </summary>
+        /// <param name="power">A number specifying a power.</param>
+        /// <returns>Returns <code>exp(power)-1</code>.</returns>
+        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;
+                }
+                );
+        }
+
+        /// <summary>
+        /// Numerically stable hypotenuse of a right angle triangle, i.e. <code>(a,b) -> sqrt(a^2 + b^2)</code>
+        /// </summary>
+        /// <param name="a">The length of side a of the triangle.</param>
+        /// <param name="b">The length of side b of the triangle.</param>
+        /// <returns>Returns <code>sqrt(a<sup>2</sup> + b<sup>2</sup>)</code> without underflow/overflow.</returns>
+        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;
+        }
+
+        /// <summary>
+        /// Numerically stable series summation
+        /// </summary>
+        /// <param name="nextSummand">provides the summands sequentially</param>
+        /// <returns>Sum</returns>
+        private static double Series(Func<double> 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;
+        }
+    }
+}

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)]