Updating Gamma distribution fitting to use Matrix's IsRelativelyEqual…

… method for determining convergence. Fixes issue accord-netGH-159.

Sources/Accord.Math/Matrix/Matrix.Common.cs

@@ -27,10 +27,13 @@ namespace Accord.Math
     using Accord.Math.Comparers;
     using System.Collections.Generic;
     using System.Collections;
+    using System.Runtime.CompilerServices;
 
     public static partial class Matrix
     {
 
+        // TODO: Use T4 templates for the equality comparisons
+
         #region Comparison
 
         /// <summary>
@@ -60,19 +63,25 @@ public static bool IsInteger(this double x, double threshold)
         ///   Compares two values for equality, considering a relative acceptance threshold.
         /// </summary>
         /// 
+#if NET45
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+#endif
         public static bool IsRelativelyEqual(this double a, double b, double threshold)
         {
-            if (Double.IsNaN(a) ^ Double.IsNaN(b))
-                return false;
+            if (a == b)
+                return true;
 
-            if (Double.IsPositiveInfinity(a) ^ Double.IsPositiveInfinity(b))
+            if (Double.IsNaN(a))
+                return Double.IsNaN(b);
+
+            if (Double.IsNaN(b)) 
                 return false;
-
-            if (Double.IsNegativeInfinity(a) ^ Double.IsNegativeInfinity(b))
+            
+            if (Double.IsPositiveInfinity(a) != Double.IsPositiveInfinity(b))
                 return false;
 
-            if (a == b)
-                return true;
+            if (Double.IsNegativeInfinity(a) != Double.IsNegativeInfinity(b))
+                return false;
 
             double delta = Math.Abs(a - b);
 

Sources/Accord.Statistics/Analysis/DistributionAnalysis.cs

@@ -140,6 +140,8 @@ public DistributionAnalysis(double[] observations)
         /// 
         public void Compute()
         {
+            bool[] fail = new bool[Distributions.Length];
+
             // Step 1. Fit all candidate distributions to the data.
             for (int i = 0; i < Distributions.Length; i++)
             {
@@ -152,6 +154,7 @@ public void Compute()
                 catch
                 {
                     // TODO: Maybe revisit the decision to swallow exceptions here.
+                    fail[i] = true;
                 }
             }
 
@@ -176,7 +179,7 @@ public void Compute()
 
                 var d = this.Distributions[i] as IUnivariateDistribution;
 
-                if (d == null)
+                if (d == null || fail[i])
                     continue;
 
                 this.DistributionNames[i] = GetName(d.GetType());

Sources/Accord.Statistics/Distributions/Univariate/Continuous/GammaDistribution.cs

@@ -439,6 +439,10 @@ public override void Fit(double[] observations, double[] weights, IFittingOption
             if (weights != null)
                 throw new ArgumentException("This distribution does not support weighted samples.");
 
+            // Method from Choi, S.C.; Wette, R. (1969) "Maximum Likelihood Estimation 
+            // of the Parameters of the Gamma Distribution and Their Bias", Technometrics,
+            // 11(4) 683–690
+
             double lnsum = 0;
             for (int i = 0; i < observations.Length; i++)
                 lnsum += Math.Log(observations[i]);
@@ -453,15 +457,17 @@ public override void Fit(double[] observations, double[] weights, IFittingOption
             // initial approximation
             double newK = (3 - s + Math.Sqrt((s - 3) * (s - 3) + 24 * s)) / (12 * s);
 
-            // Use Newton-Raphson approximation
+            // Use Newton-Raphson update
             double oldK;
 
             do
             {
                 oldK = newK;
-                newK = oldK - (Math.Log(newK) - Gamma.Digamma(newK) - s) / ((1 / newK) - Gamma.Trigamma(newK));
+                double num = Math.Log(newK) - Gamma.Digamma(newK) - s;
+                double den = (1 / newK) - Gamma.Trigamma(newK);
+                newK = oldK - num / den;
             }
-            while (Math.Abs(oldK - newK) / Math.Abs(oldK) < Double.Epsilon);
+            while (!oldK.IsRelativelyEqual(newK, 1e-10));
 
             double theta = mean / newK;
 

Unit Tests/Accord.Tests.Math/Functions/GammaTest.cs

@@ -31,22 +31,6 @@ public class GammaTest
     {
 
 
-        private TestContext testContextInstance;
-
-        public TestContext TestContext
-        {
-            get
-            {
-                return testContextInstance;
-            }
-            set
-            {
-                testContextInstance = value;
-            }
-        }
-
-
-
         [Test]
         public void FunctionTest()
         {
@@ -305,5 +289,6 @@ public void GammaUpperRTest()
             Assert.AreEqual(expected, actual, 1e-6);
             Assert.IsFalse(double.IsNaN(actual));
         }
+
     }
 }

Unit Tests/Accord.Tests.Statistics/Analysis/DistributionAnalysisTest.cs

@@ -50,6 +50,9 @@ public TestContext TestContext
 
 
         [Test]
+#if !DEBUG
+        [Timeout(2000)]
+#endif
         public void ConstructorTest()
         {
             int n = 10000;

Unit Tests/Accord.Tests.Statistics/Distributions/Univariate/Continuous/GammaDistributionTest.cs

@@ -31,22 +31,6 @@ namespace Accord.Tests.Statistics
     public class GammaDistributionTest
     {
 
-        private TestContext testContextInstance;
-
-        public TestContext TestContext
-        {
-            get
-            {
-                return testContextInstance;
-            }
-            set
-            {
-                testContextInstance = value;
-            }
-        }
-
-
-
 
         [Test]
         public void GammaDistributionConstructorTest()
@@ -235,5 +219,21 @@ public void NegativeValueTest()
 
             }
         }
+
+
+        [Test]
+        public void FitTest()
+        {
+            // Gamma Distribution Fit stalls for some arrays #159
+            // https://github.com/accord-net/framework/issues/159
+
+            double[] x = { 1275.56, 1239.44, 1237.92, 1237.22, 1237.1, 1238.41, 1238.62, 1237.05, 1237.19, 1236.51, 1264.6, 1238.19, 1237.39, 1235.79, 1236.53, 1236.8, 1238.06, 1236.5, 1235.32, 1236.44, 1236.58, 1236.3, 1237.91, 1238.6, 1238.49, 1239.21, 1238.57, 1244.63, 1236.06, 1236.4, 1237.88, 1237.56, 1236.66, 1236.59, 1236.53, 1236.32, 1238.29, 1237.79, 1237.86, 1236.42, 1236.23, 1236.37, 1237.18, 1237.63, 1245.8, 1238.04, 1238.55, 1238.39, 1236.75, 1237.07, 1250.78, 1238.6, 1238.36, 1236.58, 1236.82, 1238.4, 1257.68, 1237.78, 1236.52, 1234.9, 1237.9, 1238.58, 1238.12, 1237.89, 1236.54, 1236.55, 1238.37, 1237.29, 1237.64, 1236.8, 1237.73, 1236.71, 1238.23, 1237.84, 1236.26, 1237.58, 1238.31, 1238.4, 1237.08, 1236.61, 1235.92, 1236.41, 1237.89, 1237.98, 1246.75, 1237.92, 1237.1, 1237.97, 1238.69, 1237.05, 1236.96, 1239.44, 1238.49, 1237.88, 1236.01, 1236.57, 1236.44, 1235.76, 1237.62, 1238, 1263.14, 1237.66, 1237, 1236, 1261.96, 1238.58, 1237.77, 1237.06, 1236.31, 1238.63, 1237.23, 1236.85, 1236.23, 1236.46, 1236.9, 1237.85, 1238, 1237.02, 1236.19, 1236.05, 1235.73, 1258.3, 1235.98, 1237.76, 1246.93, 1239.1, 1237.72, 1237.67, 1236.79, 1237.61, 1238.41, 1238.29, 1238.11, 1237, 1236.52, 1236.6, 1236.31, 1237.77, 1238.58, 1237.88, 1247.35, 1236.14, 1236.83, 1236.15, 1237.93, 1238.16, 1237.34, 1236.78, 1238.66, 1237.76, 1237.19, 1236.7, 1236.04, 1236.66, 1237.86, 1238.54, 1238.05, 1238.41, 1236.94, 1240.95, 1261.01, 1237.72, 1237.91, 1238.2, 1235.68, 1236.89, 1235.12, 1271.31, 1236.97, 1270.76, 1238.52, 1238.19, 1238.6, 1237.16, 1236.72, 1236.71, 1237.14, 1238.48, 1237.95, 1237.42, 1235.86, 1236.39, 1236.13, 1236.58, 1237.95, 1237.76, 1237.39, 1238.16, 1236.31, 1236.41, 1236.12, 1238.7, 1236.48, 1237.84, 1236.38, 1237.95, 1238.48, 1236.51, 1236.56 };
+
+            var gamma = GammaDistribution.Estimate(x);
+
+            Assert.AreEqual(1238.8734170854279, gamma.Mean);
+            Assert.AreEqual(41566.439533445438, gamma.Shape);
+            Assert.AreEqual(0.029804655654680219, gamma.Scale);
+        }
     }
 }