From 011096aa8ddfd1537983c803321578e8be3a8146 Mon Sep 17 00:00:00 2001 From: Christoph Ruegg Date: Thu, 15 Feb 2018 20:43:08 +0100 Subject: [PATCH] Curve Fitting: add Fit.Exponential, Fit.Logarithm, Fit.Power --- src/Numerics.Tests/FitTests.cs | 121 +++++++++++++++++++++++++++++++++ src/Numerics/Fit.cs | 70 +++++++++++++++++++ 2 files changed, 191 insertions(+) diff --git a/src/Numerics.Tests/FitTests.cs b/src/Numerics.Tests/FitTests.cs index 0a1b7808..e4e4eb70 100644 --- a/src/Numerics.Tests/FitTests.cs +++ b/src/Numerics.Tests/FitTests.cs @@ -98,6 +98,127 @@ namespace MathNet.Numerics.UnitTests Assert.AreEqual(-0.467791, respSeq, 1e-4); } + [Test] + public void FitsToLineSameAsExcelTrendLine() + { + // X Y + // 1 0.2 + // 2 0.3 + // 4 1.3 + // 6 4.2 + // -> y = -1.078 + 0.7932*x + + var x = new[] { 1.0, 2.0, 4.0, 6.0 }; + var y = new[] { 0.2, 0.3, 1.3, 4.2 }; + + var resp = Fit.Line(x, y); + Assert.AreEqual(-1.078, resp.Item1, 1e-3); + Assert.AreEqual(0.7932, resp.Item2, 1e-3); + + var resf = Fit.LineFunc(x, y); + foreach (var z in Enumerable.Range(-3, 10)) + { + Assert.AreEqual(-1.078 + 0.7932*z, resf(z), 1e-2); + } + } + + [Test] + public void FitsToExponentialSameAsExcelTrendLine() + { + // X Y + // 1 0.2 + // 2 0.3 + // 4 1.3 + // 6 4.2 + // -> y = 0.0981*exp(0.6284*x) + + var x = new[] { 1.0, 2.0, 4.0, 6.0 }; + var y = new[] { 0.2, 0.3, 1.3, 4.2 }; + + var resp = Fit.Exponential(x, y); + Assert.AreEqual(0.0981, resp.Item1, 1e-3); + Assert.AreEqual(0.6284, resp.Item2, 1e-3); + + var resf = Fit.ExponentialFunc(x, y); + foreach (var z in Enumerable.Range(-3, 10)) + { + Assert.AreEqual(0.0981 * Math.Exp(0.6284 * z), resf(z), 1e-2); + } + } + + [Test] + public void FitsToLogarithmSameAsExcelTrendLine() + { + // X Y + // 1 0.2 + // 2 0.3 + // 4 1.3 + // 6 4.2 + // -> y = -0.4338 + 1.9981*ln(x) + + var x = new[] { 1.0, 2.0, 4.0, 6.0 }; + var y = new[] { 0.2, 0.3, 1.3, 4.2 }; + + var resp = Fit.Logarithm(x, y); + Assert.AreEqual(-0.4338, resp.Item1, 1e-3); + Assert.AreEqual(1.9981, resp.Item2, 1e-3); + + var resf = Fit.LogarithmFunc(x, y); + foreach (var z in Enumerable.Range(-3, 10)) + { + Assert.AreEqual(-0.4338 + 1.9981 * Math.Log(z), resf(z), 1e-2); + } + } + + [Test] + public void FitsToPowerSameAsExcelTrendLine() + { + // X Y + // 1 0.2 + // 2 0.3 + // 4 1.3 + // 6 4.2 + // -> y = 0.1454*x^1.7044 + + var x = new[] { 1.0, 2.0, 4.0, 6.0 }; + var y = new[] { 0.2, 0.3, 1.3, 4.2 }; + + var resp = Fit.Power(x, y); + Assert.AreEqual(0.1454, resp.Item1, 1e-3); + Assert.AreEqual(1.7044, resp.Item2, 1e-3); + + var resf = Fit.PowerFunc(x, y); + foreach (var z in Enumerable.Range(-3, 10)) + { + Assert.AreEqual(0.1454 * Math.Pow(z, 1.7044), resf(z), 1e-2); + } + } + + [Test] + public void FitsToOrder2PolynomialSameAsExcelTrendLine() + { + // X Y + // 1 0.2 + // 2 0.3 + // 4 1.3 + // 6 4.2 + // -> y = 0.7101 - 0.6675*x + 0.2077*x^2 + + var x = new[] { 1.0, 2.0, 4.0, 6.0 }; + var y = new[] { 0.2, 0.3, 1.3, 4.2 }; + + var resp = Fit.Polynomial(x, y, 2); + Assert.AreEqual(0.7101, resp[0], 1e-3); + Assert.AreEqual(-0.6675, resp[1], 1e-3); + Assert.AreEqual(0.2077, resp[2], 1e-3); + + var resf = Fit.PolynomialFunc(x, y, 2); + foreach (var z in Enumerable.Range(-3, 10)) + { + Assert.AreEqual(0.7101 - 0.6675*z + 0.2077*z*z, resf(z), 1e-2); + } + } + [Test] public void FitsToMeanOnOrder0Polynomial() { diff --git a/src/Numerics/Fit.cs b/src/Numerics/Fit.cs index 28805816..2aeb3312 100644 --- a/src/Numerics/Fit.cs +++ b/src/Numerics/Fit.cs @@ -111,6 +111,76 @@ namespace MathNet.Numerics return WeightedRegression.Weighted(x, y, w); } + /// + /// Least-Squares fitting the points (x,y) to an exponential y : x -> a*exp(r*x), + /// returning its best fitting parameters as (a, r) tuple. + /// + public static Tuple Exponential(double[] x, double[] y, DirectRegressionMethod method = DirectRegressionMethod.QR) + { + // Transformation: y_h := ln(y) ~> y_h : x -> ln(a) + r*x; + double[] y_hat = Generate.Map(y, Math.Log); + double[] p_hat = Fit.LinearCombination(x, y_hat, method, t => 1.0, t => t); + return Tuple.Create(Math.Exp(p_hat[0]), p_hat[1]); + } + + /// + /// Least-Squares fitting the points (x,y) to an exponential y : x -> a*exp(r*x), + /// returning a function y' for the best fitting line. + /// + public static Func ExponentialFunc(double[] x, double[] y, DirectRegressionMethod method = DirectRegressionMethod.QR) + { + var parameters = Exponential(x, y, method); + var a = parameters.Item1; + var r = parameters.Item2; + return z => a * Math.Exp(r * z); + } + + /// + /// Least-Squares fitting the points (x,y) to a logarithm y : x -> a + b*ln(x), + /// returning its best fitting parameters as (a, b) tuple. + /// + public static Tuple Logarithm(double[] x, double[] y, DirectRegressionMethod method = DirectRegressionMethod.QR) + { + double[] lnx = Generate.Map(x, Math.Log); + double[] p = Fit.LinearCombination(lnx, y, method, t => 1.0, t => t); + return Tuple.Create(p[0], p[1]); + } + + /// + /// Least-Squares fitting the points (x,y) to a logarithm y : x -> a + b*ln(x), + /// returning a function y' for the best fitting line. + /// + public static Func LogarithmFunc(double[] x, double[] y, DirectRegressionMethod method = DirectRegressionMethod.QR) + { + var parameters = Logarithm(x, y, method); + var a = parameters.Item1; + var b = parameters.Item2; + return z => a + b * Math.Log(z); + } + + /// + /// Least-Squares fitting the points (x,y) to a power y : x -> a*x^b, + /// returning its best fitting parameters as (a, b) tuple. + /// + public static Tuple Power(double[] x, double[] y, DirectRegressionMethod method = DirectRegressionMethod.QR) + { + // Transformation: y_h := ln(y) ~> y_h : x -> ln(a) + b*ln(x); + double[] y_hat = Generate.Map(y, Math.Log); + double[] p_hat = Fit.LinearCombination(x, y_hat, method, t => 1.0, Math.Log); + return Tuple.Create(Math.Exp(p_hat[0]), p_hat[1]); + } + /// + /// Least-Squares fitting the points (x,y) to a power y : x -> a*x^b, + /// returning a function y' for the best fitting line. + /// + public static Func PowerFunc(double[] x, double[] y, DirectRegressionMethod method = DirectRegressionMethod.QR) + { + var parameters = Power(x, y, method); + var a = parameters.Item1; + var b = parameters.Item2; + return z => a * Math.Pow(z, b); + } + /// /// 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.