From 324082c90de5637a1745735bc993bc2ee19cd56c Mon Sep 17 00:00:00 2001
From: Christoph Ruegg <git@cdrnet.ch>
Date: Sat, 29 Jun 2013 11:06:25 +0200
Subject: [PATCH] Fitting: faster algorithm for fitting to a line

---
 src/Numerics/Fit.cs | 31 ++++++++++++++++++-------------
 1 file changed, 18 insertions(+), 13 deletions(-)

diff --git a/src/Numerics/Fit.cs b/src/Numerics/Fit.cs
index fa305761..46965db5 100644
--- a/src/Numerics/Fit.cs
+++ b/src/Numerics/Fit.cs
@@ -33,7 +33,6 @@ using System.Linq;
 using MathNet.Numerics.LinearAlgebra.Double;
 using MathNet.Numerics.LinearAlgebra.Generic.Factorization;
 using MathNet.Numerics.Properties;
-using MathNet.Numerics.Statistics;
 
 namespace MathNet.Numerics
 {
@@ -53,22 +52,28 @@ namespace MathNet.Numerics
             if (x.Length != y.Length) throw new ArgumentException(Resources.ArgumentVectorsSameLength);
             if (x.Length <= 1) throw new ArgumentException(string.Format(Resources.ArrayTooSmall, 2));
 
-            var mx = ArrayStatistics.Mean(x);
-            var my = ArrayStatistics.Mean(y);
-
-            double xsum = x[0];
-            double xvariance = 0;
-            double covariance = (x[0] - mx)*(y[0] - my);
-            for (int i = 1; i < x.Length; i++)
+            // First Pass: Mean (Less robust but faster than ArrayStatistics.Mean)
+            double mx = 0.0;
+            double my = 0.0;
+            for (int i = 0; i < x.Length; i++)
             {
-                covariance += (x[i] - mx)*(y[i] - my);
+                mx += x[i];
+                my += y[i];
+            }
+            mx /= x.Length;
+            my /= y.Length;
 
-                xsum += x[i];
-                double diff = (i + 1)*x[i] - xsum;
-                xvariance += (diff*diff)/((i + 1)*i);
+            // Second Pass: Covariance/Variance
+            double covariance = 0.0;
+            double variance = 0.0;
+            for (int i = 0; i < x.Length; i++)
+            {
+                double diff = x[i] - mx;
+                covariance += diff*(y[i] - my);
+                variance += diff*diff;
             }
 
-            var b = covariance/xvariance;
+            var b = covariance/variance;
             return new[] {my - b*mx, b};
 
             // General Solution: