From dcade5cb90ab3f9406e2ee0b2f155986f3ec4768 Mon Sep 17 00:00:00 2001
From: Christoph Ruegg <git@cdrnet.ch>
Date: Wed, 26 Jun 2013 13:30:38 +0200
Subject: [PATCH] Fitting: New Least Squares Linear Curve Fitting

---
 .../FSharpPortableUnitTests.fsproj            |   3 +
 src/FSharpUnitTests/CurveFittingTests.fs      |  37 ++++
 src/FSharpUnitTests/FSharpUnitTests.fsproj    |   1 +
 src/Numerics/LeastSquares.cs                  | 110 ++++++++++++
 src/Numerics/Numerics.csproj                  |   1 +
 src/Portable/Portable.csproj                  |   3 +
 src/UnitTests/LeastSquaresTests.cs            | 160 ++++++++++++++++++
 src/UnitTests/UnitTests.csproj                |   1 +
 8 files changed, 316 insertions(+)
 create mode 100644 src/FSharpUnitTests/CurveFittingTests.fs
 create mode 100644 src/Numerics/LeastSquares.cs
 create mode 100644 src/UnitTests/LeastSquaresTests.cs

diff --git a/src/FSharpPortableUnitTests/FSharpPortableUnitTests.fsproj b/src/FSharpPortableUnitTests/FSharpPortableUnitTests.fsproj
index 6d4ad457..f106b5d3 100644
--- a/src/FSharpPortableUnitTests/FSharpPortableUnitTests.fsproj
+++ b/src/FSharpPortableUnitTests/FSharpPortableUnitTests.fsproj
@@ -91,6 +91,9 @@
     <Compile Include="..\FSharpUnitTests\PokerTests.fs">
       <Link>PokerTests.fs</Link>
     </Compile>
+    <Compile Include="..\FSharpUnitTests\CurveFittingTests.fs">
+      <Link>CurveFittingTests.fs</Link>
+    </Compile>
     <None Include="App.config" />
     <None Include="packages.config" />
   </ItemGroup>
diff --git a/src/FSharpUnitTests/CurveFittingTests.fs b/src/FSharpUnitTests/CurveFittingTests.fs
new file mode 100644
index 00000000..3e5c2a6e
--- /dev/null
+++ b/src/FSharpUnitTests/CurveFittingTests.fs
@@ -0,0 +1,37 @@
+namespace MathNet.Numerics.Tests
+
+open System
+open MathNet.Numerics
+open NUnit.Framework
+open FsUnit
+
+module CurveFittingTests =
+
+    let tofs (f:Func<_,_>) = fun a -> f.Invoke(a)
+
+    [<Test>]
+    let ``When fitting to an exact line should return exact parameters``() =
+        let f z = 4.0 - 1.5*z
+        let x = Array.append [| 1.0 .. 2.0 .. 10.0 |] [| -1.0 .. -1.0 .. -5.0  |]
+        let y = x |> Array.map f
+
+        LeastSquares.FitToLine(x,y)
+        |> should (equalWithin 1.0e-12) [| 4.0; -1.5 |]
+
+        let fres = LeastSquares.FitToLineFunc(x,y) |> tofs
+        in x |> Array.iter (fun x -> fres x |> should (equalWithin 1.0e-12) (f x))
+
+    [<Test>]
+    let ``Can fit to arbitrary linear combination``() =
+
+        // Mathematica: Fit[{{1,4.986},{2,2.347},{3,2.061},{4,-2.995},{5,-2.352},{6,-5.782}}, {1, sin(x), cos(x)}, x]
+        // -> 4.02159 sin(x) - 1.46962 cos(x) - 0.287476
+
+        let x = [| 1.0 .. 6.0 |]
+        let y = [| 4.986; 2.347; 2.061; -2.995; -2.352; -5.782 |]
+
+        LeastSquares.FitToLinearCombination(x, y, (fun z -> 1.0), (fun z -> Math.Sin(z)), (fun z -> Math.Cos(z)))
+        |> should (equalWithin 1.0e-4) [| -0.287476; 4.02159; -1.46962 |]
+
+        let fres = LeastSquares.FitToLinearCombinationFunc(x, y, (fun z -> 1.0), (fun z -> Math.Sin(z)), (fun z -> Math.Cos(z))) |> tofs
+        in x |> Array.iter (fun x -> fres x |> should (equalWithin 1.0e-4) (4.02159*Math.Sin(x) - 1.46962*Math.Cos(x) - 0.287476))
diff --git a/src/FSharpUnitTests/FSharpUnitTests.fsproj b/src/FSharpUnitTests/FSharpUnitTests.fsproj
index 75885c55..b10a7f7b 100644
--- a/src/FSharpUnitTests/FSharpUnitTests.fsproj
+++ b/src/FSharpUnitTests/FSharpUnitTests.fsproj
@@ -76,6 +76,7 @@
     <Compile Include="BigRationalTests.fs" />
     <Compile Include="RandomVariableTests.fs" />
     <Compile Include="PokerTests.fs" />
+    <Compile Include="CurveFittingTests.fs" />
     <None Include="packages.config" />
   </ItemGroup>
   <!-- To modify your build process, add your task inside one of the targets below and uncomment it. 
diff --git a/src/Numerics/LeastSquares.cs b/src/Numerics/LeastSquares.cs
new file mode 100644
index 00000000..5e7bdf56
--- /dev/null
+++ b/src/Numerics/LeastSquares.cs
@@ -0,0 +1,110 @@
+// <copyright file="LeastSquares.cs" company="Math.NET">
+// Math.NET Numerics, part of the Math.NET Project
+// http://numerics.mathdotnet.com
+// http://github.com/mathnet/mathnet-numerics
+// http://mathnetnumerics.codeplex.com
+// 
+// Copyright (c) 2009-2013 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>
+
+using System;
+using System.Linq;
+using MathNet.Numerics.LinearAlgebra.Double;
+using MathNet.Numerics.LinearAlgebra.Generic.Factorization;
+
+namespace MathNet.Numerics
+{
+    public static class LeastSquares
+    {
+        /// <summary>
+        /// Least-Squares fitting the points (x,y) to a line y : x -> a+b*x,
+        /// returning its best fitting parameters as [a, b] array.
+        /// </summary>
+        public static double[] FitToLine(double[] x, double[] y)
+        {
+            // TODO: we should use a direct algorithm instead (PERF)
+            return DenseMatrix
+                .OfColumns(x.Length, 2, new[] {DenseVector.Create(x.Length, i => 1.0), new DenseVector(x)})
+                .QR(QRMethod.Thin).Solve(new DenseVector(y))
+                .ToArray();
+        }
+
+        /// <summary>
+        /// Least-Squares fitting the points (x,y) to a line y : x -> a+b*x,
+        /// returning a function y' for the best fitting line.
+        /// </summary>
+        public static Func<double, double> FitToLineFunc(double[] x, double[] y)
+        {
+            var parameters = FitToLine(x, y);
+            double a = parameters[0], b = parameters[1];
+            return z => a + b*z;
+        }
+
+        /// <summary>
+        /// Least-Squares fitting the points (x,y) to a k-order polynomial y : x -> p0 + p1*x + p2*x^2 + ... + pk*x^k,
+        /// returning its best fitting parameters as [p0, p1, p2, ..., pk] array, compatible with Evaluate.Polynomial.
+        /// </summary>
+        public static double[] FitToPolynomial(double[] x, double[] y, int order)
+        {
+            // TODO: consider to use a specific algorithm instead
+            return DenseMatrix
+                .OfColumns(x.Length, order + 1, Enumerable.Range(0, order + 1).Select(j => DenseVector.Create(x.Length, i => Math.Pow(x[i], j))))
+                .QR(QRMethod.Thin).Solve(new DenseVector(y))
+                .ToArray();
+        }
+
+        /// <summary>
+        /// Least-Squares fitting the points (x,y) to a k-order polynomial y : x -> p0 + p1*x + p2*x^2 + ... + pk*x^k,
+        /// returning a function y' for the best fitting polynomial.
+        /// </summary>
+        public static Func<double, double> FitToPolynomialFunc(double[] x, double[] y, int order)
+        {
+            var parameters = FitToPolynomial(x, y, order);
+            return z => Evaluate.Polynomial(z, parameters);
+        }
+
+        /// <summary>
+        /// Least-Squares fitting the points (x,y) to an arbitrary linear combination y : x -> p0*f0(x) + p1*f1(x) + ... + pk*fk(x),
+        /// returning its best fitting parameters as [p0, p1, p2, ..., pk] array.
+        /// </summary>
+        public static double[] FitToLinearCombination(double[] x, double[] y, params Func<double,double>[] functions)
+        {
+            // TODO: consider to use a specific algorithm instead
+            return DenseMatrix
+                .OfColumns(x.Length, functions.Length, functions.Select(f => DenseVector.Create(x.Length, i => f(x[i]))))
+                .QR(QRMethod.Thin).Solve(new DenseVector(y))
+                .ToArray();
+        }
+
+        /// <summary>
+        /// LLeast-Squares fitting the points (x,y) to an arbitrary linear combination y : x -> p0*f0(x) + p1*f1(x) + ... + pk*fk(x),
+        /// returning a function y' for the best fitting combination.
+        /// </summary>
+        public static Func<double, double> FitToLinearCombinationFunc(double[] x, double[] y, params Func<double, double>[] functions)
+        {
+            var parameters = FitToLinearCombination(x, y, functions);
+            return z => functions.Zip(parameters, (f, p) => p*f(z)).Sum();
+        }
+    }
+}
diff --git a/src/Numerics/Numerics.csproj b/src/Numerics/Numerics.csproj
index 9f659b72..d4ecee38 100644
--- a/src/Numerics/Numerics.csproj
+++ b/src/Numerics/Numerics.csproj
@@ -105,6 +105,7 @@
     <Compile Include="Constants.cs" />
     <Compile Include="Control.cs" />
     <Compile Include="Complex32.cs" />
+    <Compile Include="LeastSquares.cs" />
     <Compile Include="Financial\AbsoluteReturnMeasures.cs" />
     <Compile Include="Financial\AbsoluteRiskMeasures.cs" />
     <Compile Include="LinearAlgebra\Generic\Matrix.BCL.cs" />
diff --git a/src/Portable/Portable.csproj b/src/Portable/Portable.csproj
index b7c8bed7..4cf36a4f 100644
--- a/src/Portable/Portable.csproj
+++ b/src/Portable/Portable.csproj
@@ -282,6 +282,9 @@
     <Compile Include="..\Numerics\IPrecisionSupport.cs">
       <Link>IPrecisionSupport.cs</Link>
     </Compile>
+    <Compile Include="..\Numerics\LeastSquares.cs">
+      <Link>LeastSquares.cs</Link>
+    </Compile>
     <Compile Include="..\Numerics\LinearAlgebra\Complex32\DenseMatrix.cs">
       <Link>LinearAlgebra\Complex32\DenseMatrix.cs</Link>
     </Compile>
diff --git a/src/UnitTests/LeastSquaresTests.cs b/src/UnitTests/LeastSquaresTests.cs
new file mode 100644
index 00000000..176cfd83
--- /dev/null
+++ b/src/UnitTests/LeastSquaresTests.cs
@@ -0,0 +1,160 @@
+// <copyright file="LeastSquaresTests.cs" company="Math.NET">
+// Math.NET Numerics, part of the Math.NET Project
+// http://numerics.mathdotnet.com
+// http://github.com/mathnet/mathnet-numerics
+// http://mathnetnumerics.codeplex.com
+// 
+// Copyright (c) 2009-2013 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>
+
+using System;
+using System.Linq;
+using MathNet.Numerics.Statistics;
+using NUnit.Framework;
+
+namespace MathNet.Numerics.UnitTests
+{
+    [TestFixture]
+    public class LeastSquaresTests
+    {
+        [Test]
+        public void FitsToExactLineWhenPointsAreOnLine()
+        {
+            var x = new[] {30.0, 40.0, 50.0, 12.0, -3.4, 100.5};
+            var y = x.Select(z => 4.0 - 1.5*z).ToArray();
+
+            var resp = LeastSquares.FitToLine(x, y);
+            Assert.AreEqual(2, resp.Length);
+            Assert.AreEqual(4.0, resp[0], 1e-12);
+            Assert.AreEqual(-1.5, resp[1], 1e-12);
+
+            var resf = LeastSquares.FitToLineFunc(x, y);
+            foreach (var z in Enumerable.Range(-3, 10))
+            {
+                Assert.AreEqual(4.0 - 1.5*z, resf(z), 1e-12);
+            }
+        }
+
+        [Test]
+        public void FitsToBestLine()
+        {
+            // Mathematica: Fit[{{1,4.986},{2,2.347},{3,2.061},{4,-2.995},{5,-2.352},{6,-5.782}}, {1, x}, x]
+            // -> 7.01013 - 2.08551*x
+
+            var x = Enumerable.Range(1, 6).Select(Convert.ToDouble).ToArray();
+            var y = new[] {4.986, 2.347, 2.061, -2.995, -2.352, -5.782};
+
+            var resp = LeastSquares.FitToLine(x, y);
+            Assert.AreEqual(2, resp.Length);
+            Assert.AreEqual(7.01013, resp[0], 1e-4);
+            Assert.AreEqual(-2.08551, resp[1], 1e-4);
+
+            var resf = LeastSquares.FitToLineFunc(x, y);
+            foreach (var z in Enumerable.Range(-3, 10))
+            {
+                Assert.AreEqual(7.01013 - 2.08551 * z, resf(z), 1e-4);
+            }
+        }
+
+        [Test]
+        public void FitsToMeanOnOrder0Polynomial()
+        {
+            // Mathematica: Fit[{{1,4.986},{2,2.347},{3,2.061},{4,-2.995},{5,-2.352},{6,-5.782}}, {1}, x]
+            // -> -0.289167
+
+            var x = Enumerable.Range(1, 6).Select(Convert.ToDouble).ToArray();
+            var y = new[] { 4.986, 2.347, 2.061, -2.995, -2.352, -5.782 };
+
+            var resp = LeastSquares.FitToPolynomial(x, y, 0);
+            Assert.AreEqual(1, resp.Length);
+            Assert.AreEqual(-0.289167, resp[0], 1e-4);
+            Assert.AreEqual(y.Mean(), resp[0], 1e-4);
+        }
+
+        [Test]
+        public void FitsToLineOnOrder1Polynomial()
+        {
+            // Mathematica: Fit[{{1,4.986},{2,2.347},{3,2.061},{4,-2.995},{5,-2.352},{6,-5.782}}, {1, x}, x]
+            // -> 7.01013 - 2.08551 x
+
+            var x = Enumerable.Range(1, 6).Select(Convert.ToDouble).ToArray();
+            var y = new[] { 4.986, 2.347, 2.061, -2.995, -2.352, -5.782 };
+
+            var resp = LeastSquares.FitToPolynomial(x, y, 1);
+            Assert.AreEqual(2, resp.Length);
+            Assert.AreEqual(7.01013, resp[0], 1e-4);
+            Assert.AreEqual(-2.08551, resp[1], 1e-4);
+
+            var resf = LeastSquares.FitToPolynomialFunc(x, y, 1);
+            foreach (var z in Enumerable.Range(-3, 10))
+            {
+                Assert.AreEqual(7.01013 - 2.08551 * z, resf(z), 1e-4);
+            }
+        }
+
+        [Test]
+        public void FitsToOrder2Polynomial()
+        {
+            // Mathematica: Fit[{{1,4.986},{2,2.347},{3,2.061},{4,-2.995},{5,-2.352},{6,-5.782}}, {1, x, x^2}, x]
+            // -> 6.9703 - 2.05564 x - 0.00426786 x^2
+
+            var x = Enumerable.Range(1, 6).Select(Convert.ToDouble).ToArray();
+            var y = new[] { 4.986, 2.347, 2.061, -2.995, -2.352, -5.782 };
+
+            var resp = LeastSquares.FitToPolynomial(x, y, 2);
+            Assert.AreEqual(3, resp.Length);
+            Assert.AreEqual(6.9703, resp[0], 1e-4);
+            Assert.AreEqual(-2.05564, resp[1], 1e-4);
+            Assert.AreEqual(-0.00426786, resp[2], 1e-6);
+
+            var resf = LeastSquares.FitToPolynomialFunc(x, y, 2);
+            foreach (var z in Enumerable.Range(-3, 10))
+            {
+                Assert.AreEqual(Evaluate.Polynomial(z, resp), resf(z), 1e-4);
+            }
+        }
+
+        [Test]
+        public void FitsToTrigonometricLinearCombination()
+        {
+            // Mathematica: Fit[{{1,4.986},{2,2.347},{3,2.061},{4,-2.995},{5,-2.352},{6,-5.782}}, {1, sin(x), cos(x)}, x]
+            // -> 4.02159 sin(x) - 1.46962 cos(x) - 0.287476
+
+            var x = Enumerable.Range(1, 6).Select(Convert.ToDouble).ToArray();
+            var y = new[] { 4.986, 2.347, 2.061, -2.995, -2.352, -5.782 };
+
+            var resp = LeastSquares.FitToLinearCombination(x, y, z => 1.0, Math.Sin, Math.Cos);
+            Assert.AreEqual(3, resp.Length);
+            Assert.AreEqual(-0.287476, resp[0], 1e-4);
+            Assert.AreEqual(4.02159, resp[1], 1e-4);
+            Assert.AreEqual(-1.46962, resp[2], 1e-4);
+
+            var resf = LeastSquares.FitToLinearCombinationFunc(x, y, z => 1.0, Math.Sin, Math.Cos);
+            foreach (var z in Enumerable.Range(-3, 10))
+            {
+                Assert.AreEqual(4.02159*Math.Sin(z) - 1.46962*Math.Cos(z) - 0.287476, resf(z), 1e-4);
+            }
+        }
+    }
+}
diff --git a/src/UnitTests/UnitTests.csproj b/src/UnitTests/UnitTests.csproj
index 68526483..d13c9564 100644
--- a/src/UnitTests/UnitTests.csproj
+++ b/src/UnitTests/UnitTests.csproj
@@ -144,6 +144,7 @@
     <Compile Include="InterpolationTests\EquidistantPolynomialTest.cs" />
     <Compile Include="InterpolationTests\NevillePolynomialTest.cs" />
     <Compile Include="InterpolationTests\LinearInterpolationCase.cs" />
+    <Compile Include="LeastSquaresTests.cs" />
     <Compile Include="LinearAlgebraProviderTests\Complex32\LinearAlgebraProviderTests.cs" />
     <Compile Include="LinearAlgebraProviderTests\Complex\LinearAlgebraProviderTests.cs" />
     <Compile Include="LinearAlgebraProviderTests\Double\LinearAlgebraProviderTests.cs">