diff --git a/src/Numerics/GoodnessOfFit.cs b/src/Numerics/GoodnessOfFit.cs
index 5972384f..29ad6cca 100644
--- a/src/Numerics/GoodnessOfFit.cs
+++ b/src/Numerics/GoodnessOfFit.cs
@@ -3,6 +3,8 @@
 // http://numerics.mathdotnet.com
 // http://github.com/mathnet/mathnet-numerics
 //
+// Copyright (c) 2009-2018 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
@@ -61,7 +63,7 @@ namespace MathNet.Numerics
 
         /// <summary>
         /// Calculates the Standard Error of the regression, given a sequence of
-        /// modeled/predicted values, and a sequence of actual/observed values 
+        /// modeled/predicted values, and a sequence of actual/observed values
         /// </summary>
         /// <param name="modelledValues">The modelled/predicted values</param>
         /// <param name="observedValues">The observed/actual values</param>
@@ -77,7 +79,7 @@ namespace MathNet.Numerics
         /// </summary>
         /// <param name="modelledValues">The modelled/predicted values</param>
         /// <param name="observedValues">The observed/actual values</param>
-        /// <param name="degreesOfFreedom">The degrees of freedom by which the 
+        /// <param name="degreesOfFreedom">The degrees of freedom by which the
         /// number of samples is reduced for performing the Standard Error calculation</param>
         /// <returns>The Standard Error of the regression</returns>
         public static double StandardError(IEnumerable<double> modelledValues, IEnumerable<double> observedValues, int degreesOfFreedom)
@@ -107,5 +109,65 @@ namespace MathNet.Numerics
                 return Math.Sqrt(accumulator / (n - degreesOfFreedom));
             }
         }
+
+        /// <summary>
+        /// Calculates the R-Squared value, also known as coefficient of determination,
+        /// given some modelled and observed values.
+        /// </summary>
+        /// <param name="modelledValues">The values expected from the model.</param>
+        /// <param name="observedValues">The actual values obtained.</param>
+        /// <returns>Coefficient of determination.</returns>
+        public static double CoefficientOfDetermination(IEnumerable<double> modelledValues, IEnumerable<double> observedValues)
+        {
+            var y = observedValues;
+            var f = modelledValues;
+            int n = 0;
+
+            double meanY = 0;
+            double ssTot = 0;
+            double ssRes = 0;
+
+            using (IEnumerator<double> ieY = y.GetEnumerator())
+            using (IEnumerator<double> ieF = f.GetEnumerator())
+            {
+                while (ieY.MoveNext())
+                {
+                    if (!ieF.MoveNext())
+                    {
+                        throw new ArgumentOutOfRangeException("modelledValues", Resources.ArgumentArraysSameLength);
+                    }
+
+                    double currentY = ieY.Current;
+                    double currentF = ieF.Current;
+
+                    // If a large constant C is added to every y value,
+                    // then each new y have an error of about C*eps,
+                    // thus each new deltaY will change by about C*eps (compared to the old deltaY),
+                    // and thus ssTot will change by only C*eps*deltaY on each step
+                    // thus C*eps*deltaY*n in total.
+                    // (This error cannot be eliminated by a Kahan algorithm,
+                    // because it is introduced when C is added to the old Y value).
+                    //
+                    // This is better than summing the square of y values
+                    // and then substracting the correct multiple of the square of the sum of y values,
+                    // in this latter case ssTot will change by eps*n*(C^2+2*C*meanY) in total.
+                    double deltaY = currentY - meanY;
+                    double scaleDeltaY = deltaY / ++n;
+
+                    meanY += scaleDeltaY;
+                    ssTot += scaleDeltaY* deltaY* (n - 1);
+
+                    // This calculation is as safe as ssTot
+                    // in the case when a constant is added to both y and f.
+                    ssRes += (currentY - currentF)* (currentY-currentF);
+                }
+
+                if (ieF.MoveNext())
+                {
+                    throw new ArgumentOutOfRangeException("observedValues", Resources.ArgumentArraysSameLength);
+                }
+            }
+            return 1 - ssRes/ssTot;
+        }
     }
 }