--- title: "Fit" language: "en" type: "Symbol" summary: "Fit[data, {f1, ..., fn}, {x, y, ...}] finds a fit a1\\[InvisibleTimes]f1 + ... + an\\[InvisibleTimes]fn to a list of data for functions f1, ..., fn of variables {x, y, ...}. Fit[{m, v}] finds a fit vector a that minimizes || m . a - v || for a design matrix m. Fit[..., prop] specifies what fit property prop should be returned." keywords: - approximate formulas - approximation of functions - chi squared - curves - data - fits - fitting of data - least-squares fits - linear fits - models - robust fitting - huber penalty - regularized fitting - LASSO - total variation regularization - l1 regularization - l2 regularization - Tikhonov regularization - ridge regression - basis pursuit - COMFIT - CURVEFIT - GAUSS2DFIT - GAUSSFIT - LADFIT - LINFIT - LMFIT - POLY_FIT - REGRESS - SFIT - CurveFitting - linfit - polyfit canonical_url: "https://reference.wolfram.com/language/ref/Fit.html" source: "Wolfram Language Documentation" related_guides: - title: "Curve Fitting & Approximate Functions" link: "https://reference.wolfram.com/language/guide/CurveFittingAndApproximateFunctions.en.md" - title: "Data Transforms and Smoothing" link: "https://reference.wolfram.com/language/guide/DataTransformsAndSmoothing.en.md" - title: "Matrix-Based Minimization" link: "https://reference.wolfram.com/language/guide/MatrixBasedMinimization.en.md" - title: "Numerical Data" link: "https://reference.wolfram.com/language/guide/NumericalData.en.md" related_functions: - title: "FindFit" link: "https://reference.wolfram.com/language/ref/FindFit.en.md" - title: "LinearModelFit" link: "https://reference.wolfram.com/language/ref/LinearModelFit.en.md" - title: "NonlinearModelFit" link: "https://reference.wolfram.com/language/ref/NonlinearModelFit.en.md" - title: "DesignMatrix" link: "https://reference.wolfram.com/language/ref/DesignMatrix.en.md" - title: "LeastSquares" link: "https://reference.wolfram.com/language/ref/LeastSquares.en.md" - title: "PseudoInverse" link: "https://reference.wolfram.com/language/ref/PseudoInverse.en.md" - title: "Interpolation" link: "https://reference.wolfram.com/language/ref/Interpolation.en.md" - title: "InterpolatingPolynomial" link: "https://reference.wolfram.com/language/ref/InterpolatingPolynomial.en.md" - title: "Solve" link: "https://reference.wolfram.com/language/ref/Solve.en.md" - title: "ListPlot" link: "https://reference.wolfram.com/language/ref/ListPlot.en.md" - title: "ConvexOptimization" link: "https://reference.wolfram.com/language/ref/ConvexOptimization.en.md" related_tutorials: - title: "Manipulating Numerical Data" link: "https://reference.wolfram.com/language/tutorial/NumericalOperationsOnData.en.md#25465" - title: "Curve Fitting" link: "https://reference.wolfram.com/language/tutorial/NumericalOperationsOnData.en.md#19037" - title: "Implementation notes: Numerical and Related Functions" link: "https://reference.wolfram.com/language/tutorial/SomeNotesOnInternalImplementation.en.md#20880" --- # Fit Fit[data, {f1, …, fn}, {x, y, …}] finds a fit a1⁢f1 + … + an⁢fn to a list of data for functions f1, …, fn of variables {x, y, …}. Fit[{m, v}] finds a fit vector a that minimizes $\|m.a-v\|$ for a design matrix m. Fit[…, "prop"] specifies what fit property prop should be returned. ## Details and Options * ``Fit`` is also known as linear regression or least squares fit. With regularization, it is also known as LASSO and ridge regression. * ``Fit`` is typically used for fitting combinations of functions to data, including polynomials and exponentials. It provides one of the simplest ways to get a model from data. [image] * The best fit minimizes the sum of squares $\left(a_1f_1\left(x_1,y_1,\ldots \right)+\ldots +a_nf_n\left(x_1,y_1,\ldots \right)-v_1\right){}^2+\ldots +\left(a_1f_1\left(x_k,y_k,\ldots \right)+\ldots +a_nf_n\left(x_k,y_k,\ldots \right)-v_k\right){}^2$. * The ``data`` can have the following forms: | | | | ----------------------- | ----------------------------------------------------------- | | {v1, …, vn} | equivalent to {{1, v1}, …, {n, vn}} | | {{x1, v1}, …, {xn, vn}} | univariate data with values vi at coordinates xi | | {{x1, y1, v1}, …} | bivariate data with values vi and coordinates {xi, yi} | | {{x1, y1, …, v1}, …} | multivariate data with values vi at coordinates {xi, yi, …} | * The design matrix ``m`` has elements that come from evaluating the functions at the coordinates, $m_{ij}=f_i\left(x_j,y_j,\ldots \right)$. In matrix notation, the best fit minimizes the norm $\| m.a-v\|$ where $a=\left\{a_1,\ldots ,a_n\right\}$ and $v=\left\{v_1,\ldots ,v_k\right\}$. * The functions ``fi`` should depend only on the variables ``{x, y, …}``. * The possible fit properties ``"prop"`` include: "BasisFunctions" funs the basis functions "BestFit" Subscript[a, 1]\[InvisibleTimes]Subscript[f, 1]+\[Ellipsis]+Subscript[a, n]\[InvisibleTimes]Subscript[f, n] the best fit linear combination of basis functions "BestFitParameters" {Subscript[a, 1],\[Ellipsis],Subscript[a, n]} the vector a that gives the best fit "Coordinates" {{Subscript[x, 1],Subscript[y, 1],\[Ellipsis]},\[Ellipsis]} the coordinates of vars in data "Data" data the data "DesignMatrix" m (Subscript[f, 1](Subscript[x, 1],Subscript[y, 1],\[Ellipsis]) \[Ellipsis] Subscript[f, n](Subscript[x, 1],Subscript[y, 1],\[Ellipsis]) \[VerticalEllipsis] \[DescendingEllipsis] \[VerticalEllipsis] Subscript[f, 1](Subscript[x, k],Subscript[y, k],\[Ellipsis]) \[Ellipsis] Subscript[f, n](Subscript[x, k],Subscript[y, k],\[Ellipsis]) ) "FitResiduals" m.a-v the differences between the model and the fit at the coordinates "Function" Function[{x,y,\[Ellipsis]},Subscript[a, 1] Subscript[f, 1]+\[Ellipsis]+Subscript[a, n] Subscript[f, n]] best fit pure function "PredictedResponse" m.a fitted values for the data coordinates "Response" {Subscript[v, 1],\[Ellipsis],Subscript[v, k]} the response vector v from the input data {"Subscript[prop, 1]","Subscript[prop, 2]",\[Ellipsis]} several fit properties * ``Fit`` takes the following options: | | | | | ------------------ | --------- | ----------------------------------------------------------------------- | | NormFunction | Norm | measure of the deviations to minimize | | FitRegularization | None | regularization for the fit parameters $a$ | | WorkingPrecision | Automatic | the precision to use | * With ``NormFunction -> normf`` and ``FitRegularization -> rfun``, ``Fit`` finds the coefficient vector ``a`` that minimizes ``normf[{a.f(x1, y1, …) - v1, …, a.f(xk, yk, …) - vk}] + rfun[a]``. * The setting for ``NormFunction`` can be given in the following forms: | | | | ---------------------------- | -------------------------------------------------------------------------- | | normf | a function normf that is applied to the deviations | | {"Penalty", pf} | sum of the penalty function pf applied to each component of the deviations | | {"HuberPenalty", α} | sum of Huber penalty function for each component | | {"DeadzoneLinearPenalty", α} | sum of deadzone linear penalty function for each component | * The setting for ``FitRegularization`` may be given in the following forms: | | | | --- | --- | | None | no regularization | | rfun | regularize with rfun[a] | | {"Tikhonov", λ} | regularize with $\lambda \left\\|a\\|^2\right.$ | | {"LASSO", λ} | regularize with $\lambda \left\\|a\\|_1\right.$ | | {"Variation", λ} | regularize with $\lambda \left\\|\fbox{$\text{Differences}$}[a]\\|^2\right.$ | | {"TotalVariation", λ} | regularize with $\lambda \left\\|\fbox{$\text{Differences}$}[a]\\|_1\right.$ | | {"Curvature", λ} | regularize with $\lambda \left\\|\fbox{$\text{Differences}$}[a,2]\\|^2\right.$ | | {r1, r2, …} | regularize with the sum of terms from r1, … | * With ``WorkingPrecision -> Automatic``, exact numbers given as input to ``Fit`` are converted to approximate numbers with machine precision. --- ## Examples (24) ### Basic Examples (2) Here is some data: ```wl In[1]:= data = {{0, 1}, {1, 0}, {3, 2}, {5, 4}}; ``` Find the line that best fits the data: ```wl In[2]:= line = Fit[data, {1, x}, x] Out[2]= 0.186441  + 0.694915 x In[3]:= DesignMatrix[data, {1, x}, x] Out[3]= {{1, 0}, {1, 1}, {1, 3}, {1, 5}} ``` Find the quadratic that best fits the data: ```wl In[4]:= parabola = Fit[data, {1, x, x ^ 2}, x] Out[4]= 0.678392  - 0.266332 x + 0.190955 x^2 ``` Show the data with the two curves: ```wl In[5]:= Show[ListPlot[data, PlotStyle -> Red], Plot[{line, parabola}, {x, 0, 5}]] Out[5]= [image] ``` --- Find the best fit parameters given a design matrix and response vector: ```wl In[1]:= m = N[HilbertMatrix[4]]; v = Range[4]; a = Fit[{m, v}] Out[1]= {-64., 900., -2520., 1820.} ``` ### Scope (2) Here is some data defined with exact values: ```wl In[1]:= data = {{-Pi, 4}, {-Pi / 2, 0}, {0, 1}, {Pi / 2, -1}, {Pi, -4}}; ``` Fit the data to a linear combination of sine functions using machine arithmetic: ```wl In[2]:= Fit[data, {Sin[x / 2], Sin[x], Sin[2 x]}, x] Out[2]= 0.  - 4. Sin[(x/2)] + 2.32843 Sin[x] ``` Fit the data using 24-digit precision arithmetic: ```wl In[3]:= Fit[N[data, 24], {Sin[x / 2], Sin[x], Sin[2 x]}, x] Out[3]= -4.000000000000000000000 Sin[(x/2)] + 2.328427124746190097603 Sin[x] ``` Show the data with the curve: ```wl In[4]:= Show[ListPlot[data], Plot[%, {x, -Pi, Pi}], PlotRange -> All] Out[4]= [image] ``` --- Here is some data in three dimensions: ```wl In[1]:= data = {{0, 0, 0}, {1, 0, 1}, {0, 1, 2}, {1, 1, 0}, {1 / 2, 1 / 2, 1}}; ``` Find the plane that best fits the data: ```wl In[2]:= plane = Fit[data, {1, x, y}, {x, y}] Out[2]= 0.8  - 0.5 x + 0.5 y ``` Show the plane with the data points: ```wl In[3]:= Show[Plot3D[plane, {x, 0, 1}, {y, 0, 1}, PlotStyle -> Opacity[.5], PlotRange -> {0, 2}], Graphics3D[{Red, PointSize[0.05], Map[Point, data]}]] Out[3]= [image] ``` Find the quadratic that best fits the data: ```wl In[4]:= quad = Fit[data, {1, x, y, x ^ 2, x y, y ^ 2}, {x, y}] Out[4]= 3.856346010722665`*^-16 + 1.76087 x - 0.76087 x^2 + 2.23913 y - 3. x y - 0.23913 y^2 ``` The quadratic actually interpolates the data: ```wl In[5]:= Show[Plot3D[quad, {x, 0, 1}, {y, 0, 1}, PlotStyle -> Opacity[.5], PlotRange -> {0, 2}], Graphics3D[{Red, PointSize[0.05], Map[Point, data]}]] Out[5]= [image] ``` ### Generalizations & Extensions (1) Here is a list of values: ```wl In[1]:= v = {5, 4, 3, 2, 1, 0, 1, 2, 3, 4, 5}; ``` Fit to a quadratic. When coordinates are not given, the values are assumed to be paired up with 1, 2, …: ```wl In[2]:= quad = Fit[v, {1, x, x ^ 2}, x] Out[2]= 7.27273  - 2.0979 x + 0.174825 x^2 ``` Fit to a quartic: ```wl In[3]:= quart = Fit[v, {1, x, x ^ 2, x ^ 3, x ^ 4}, x] Out[3]= 4.54545  + 1.18881 x - 0.938228 x^2 + 0.13986 x^3 - 0.00582751 x^4 ``` Show the data with the curve: ```wl In[4]:= Show[ListPlot[v, PlotStyle -> Black], Plot[{quad, quart}, {x, 1, Length[v]}]] Out[4]= [image] ``` ### Options (6) #### FitRegularization (2) Tikhonov regularization controls the size of the best fit parameters: ```wl In[1]:= data = {{0., 0.}, {0.001, 1}, {0.01, 1}}; In[2]:= Fit[data, {1, x, x ^ 2}, x, FitRegularization -> {"Tikhonov", 1}] Out[2]= 0.499985  + 0.00549961 x + 0.000050496 x^2 ``` Without the regularization, the coefficients are quite large: ```wl In[3]:= Fit[data, {1, x, x ^ 2}, x] Out[3]= 1.832101568702095`*^-15 + 1100. x - 100000. x^2 ``` This can also be referred to by ``"L2"`` or ``"RidgeRegression"`` : ```wl In[4]:= Fit[data, {1, x, x ^ 2}, x, FitRegularization -> {"L2", 1}] Out[4]= 0.499985  + 0.00549961 x + 0.000050496 x^2 In[5]:= Fit[data, {1, x, x ^ 2}, x, FitRegularization -> {"RidgeRegression", 1}] Out[5]= 0.499985  + 0.00549961 x + 0.000050496 x^2 ``` --- LASSO (least absolute shrinkage and selection operator) regularization selects the most significant basis functions and gives some control of the size of the best fit parameters: ```wl In[1]:= data = {{0., 0.}, {0.001, 1}, {0.01, 1}, {.1, 0}, {1, 1}}; In[2]:= Fit[data, {1, x, x ^ 2, x ^ 3, x ^ 4}, x, FitRegularization -> {"LASSO", 1}] Out[2]= 0.499939  + 0.000229385 x^4 ``` Without the regularization, the coefficients are quite large: ```wl In[3]:= Fit[data, {1, x, x ^ 2, x ^ 3, x ^ 4}, x] Out[3]= 1.3291973418163954`*^-10 + 1110.99 x - 112096. x^2 + 1.1097392442134216*^6 x^3 - 998753. x^4 ``` This can also be referred to by ``"L1"`` : ```wl In[4]:= Fit[data, {1, x, x ^ 2, x ^ 3, x ^ 4}, x, FitRegularization -> {"L1", 1}] Out[4]= 0.499939  + 0.000229385 x^4 ``` #### NormFunction (3) Measure the "best" fit using the 1-norm: ```wl In[1]:= data = {{0, 1}, {1, 0}, {3, 2}, {5, 4}}; In[2]:= l1 = Fit[data, {1, x}, x, NormFunction -> Function[Norm[#, 1]]] Out[2]= -1. + 1. x In[3]:= Show[ListPlot[data, PlotStyle -> PointSize[0.025]], Plot[l1, {x, 0, 5}], PlotRange -> All] Out[3]= [image] ``` Compare with the least squares (gold) and max norm (green) fits: ```wl In[4]:= l2 = Fit[data, {1, x}, x]; l∞ = Fit[data, {1, x}, x, NormFunction -> Function[Norm[#, ∞]]]; Show[ListPlot[data, PlotStyle -> PointSize[0.025]], Plot[{l1, l2, l∞}, {x, 0, 5}], PlotRange -> All] Out[4]= [image] ``` Manipulate the norm ``p`` interactively: ```wl In[5]:= Manipulate[Plot[Evaluate[With[{q = 1 / pr}, Fit[data, {1, x}, x, NormFunction -> Function[Norm[#, q]]]]], {x, 0, 5}, PlotRange -> {-1, 4}, Epilog -> {PointSize[.025], Point[data], Text["p = " <> ToString[1 / pr], {1, 3}]}], {{pr, 1 / 2, "1/p"}, 0.0001, .9999}, SaveDefinitions -> True] Out[5]= DynamicModule[«9»] ``` --- Use the Huber penalty function to balance the influence of outliers with least squares: ```wl In[1]:= data = {{0, 1}, {1, 0}, {3, 2}, {5, 4}}; In[2]:= hfit = Fit[data, {1, x}, x, NormFunction -> {"HuberPenalty", .1}] Out[2]= -0.8 + 0.95 x In[3]:= Show[ListPlot[data, PlotStyle -> PointSize[0.025]], Plot[hfit, {x, 0, 5}], PlotRange -> All] Out[3]= [image] ``` Manipulate the parameter $\alpha$ interactively: ```wl In[4]:= Manipulate[Plot[Evaluate[Fit[data, {1, x}, x, NormFunction -> {"HuberPenalty", α}]], {x, 0, 5}, PlotRange -> {-1, 4}, Epilog -> {PointSize[.025], Point[data], Text["α = " <> ToString[α], {1, 3}]}], {{α, .1}, 0.0001, 2}, SaveDefinitions -> True] Out[4]= DynamicModule[«9»] ``` --- Compute a weighted least squares fit: ```wl In[1]:= data = {{0, 1}, {1, 0}, {3, 2}, {5, 4}}; fit = Fit[data, {1, x}, x, NormFunction -> Function[Norm[# Range[4]]]]; Show[ListPlot[data, PlotStyle -> PointSize[0.025]], Plot[fit, {x, 0, 5}], PlotRange -> All] Out[1]= [image] ``` #### WorkingPrecision (1) Compute the fit using 24-digit precision: ```wl In[1]:= data = {{0, 1}, {1, 0}, {3, 2}, {5, 4}}; In[2]:= Fit[data, {1, x, x ^ 2}, x, WorkingPrecision -> 24] Out[2]= 0.6783919597989949748744 - 0.2663316582914572864322 x + 0.19095477386934673366834 x^2 ``` Use exact computations: ```wl In[3]:= Fit[data, {1, x, x ^ 2}, x, WorkingPrecision -> ∞] Out[3]= (135/199) - (53 x/199) + (38 x^2/199) ``` ### Applications (6) Compute a high-degree polynomial fit to data using a Chebyshev basis: ```wl In[1]:= coord = RandomReal[{-1, 1}, 100]; data = Transpose[{coord, Sinc[10 coord]}]; In[2]:= degree = 47; basis = Table[ChebyshevT[k, x], {k, 0, degree}]; fit = Fit[data, basis, x]; In[3]:= Plot[fit, {x, -1, 1}, Epilog -> Point[data]] Out[3]= [image] ``` --- Use regularization to stabilize the numerical solution $x$ to $a.x=b$ where $b$ is perturbed: ```wl In[1]:= a = N[HilbertMatrix[47]]; xexact = Table[t * (1 - t), {t, 0., 1., 1 / 46}]; b = a.xexact + 10^ - 8 RandomReal[{-1, 1}, 47]; ``` The solution found by ``LinearSolve`` has very large entries: ```wl In[2]:= ListPlot[LinearSolve[a, b]] ``` LinearSolve::luc: Result for LinearSolve of badly conditioned matrix {{1.,0.5,0.333333,0.25,0.2,0.166667,0.142857,0.125,0.111111,0.1,0.0909091,0.0833333,<<24>>,0.027027,0.0263158,0.025641,0.025,0.0243902,0.0238095,0.0232558,0.0227273,0.0222222,0.0217391,0.0212766},<<45>>,{<<1>>}} may contain significant numerical errors. ```wl Out[2]= [image] ``` Regularized solutions stay much closer to the solution of the unperturbed problem: ```wl In[3]:= xreg = Table[Fit[{a, b}, FitRegularization -> {"Tikhonov", 10 ^ -ll}], {ll, 0, 3, 2}]; In[4]:= ListPlot[Prepend[xreg, xexact], PlotStyle -> PointSize[0.025]] Out[4]= [image] ``` --- Use variation regularization to find a smooth approximation for an input signal: ```wl In[1]:= n = 201; m = 201; times = Range[m]; ``` The output target is a step function: ```wl In[2]:= m1 = Round[m / 5]; m2 = Round[m / 4]; m3 = Round[m / 4];output = ConstantArray[0, m]; output[[m1 + Range[m2]]] = 1; output[[m1 + m2 + Range[m3]]] = -1; op = ListLinePlot[output] Out[2]= [image] ``` The input is determined from the output via convolution with $h=\frac{1}{9}\left(1-\left.\frac{2}{5}\cos (2 t)\right)\left(\frac{9}{10}\right)^t\right.$ : ```wl In[3]:= h = (1/9)(1 - .4Cos[2 times]) * .9 ^ times;hmat = ToeplitzMatrix[h, Prepend[ConstantArray[0, n - 1], h[[1]]]]; ``` Without regularization, the predicated response matches the signal very closely, but the computed input has a lot of oscillation: ```wl In[4]:= {input, predicted} = Fit[{hmat, output}, {"BestFitParameters", "PredictedResponse"}]; Row[{ListLinePlot[input], ListLinePlot[{predicted, output}]}] Out[4]= [image] ``` With regularization of the variation, a much smoother approximation can be found: ```wl In[5]:= {input, predicted} = Fit[{hmat, output}, {"BestFitParameters", "PredictedResponse"}, FitRegularization -> {"Variation", .3}]; Row[{ListLinePlot[input], ListLinePlot[{predicted, output}]}] Out[5]= [image] ``` Regularization of the input size may also be included: ```wl In[6]:= {input, predicted} = Fit[{hmat, output}, {"BestFitParameters", "PredictedResponse"}, FitRegularization -> {{"Variation", .3}, {"Tikhonov", .01}}]; Row[{ListLinePlot[input], ListLinePlot[{predicted, output}]}] Out[6]= [image] ``` --- Smooth a corrupted signal using variation regularization: ```wl In[1]:= n = 4000; times = N[Range[0, n - 1]]; original = .5 Sin[2 Pi times / n] * Sin[.01 times]; corrupted = original + RandomReal[{-.05, .05}, n]; ListLinePlot[corrupted] Out[1]= [image] ``` Plot the tradeoff between the norm of the residuals and the norm of the variation: ```wl In[2]:= id = IdentityMatrix[n, SparseArray]; In[3]:= tradeoff = Table[{smoothed, residual} = Fit[{id, corrupted}, {"BestFitParameters", "FitResiduals"}, FitRegularization -> {"Variation", 2 ^ logλ}];Tooltip[{Norm[Differences[smoothed]], Norm[residual]}, 2. ^ logλ], {logλ, -10, 20}]; In[4]:= ListPlot[tradeoff, Joined -> True, Mesh -> All] Out[4]= [image] ``` Choose a value of $\lambda$ near where the curve bends sharply: ```wl In[5]:= smoothed = Fit[{id, corrupted}, "BestFitParameters", FitRegularization -> {"Variation", 100}]; In[6]:= ListLinePlot[{corrupted, smoothed}] Out[6]= [image] ``` --- Use total variation regularization to smooth a corrupted signal with jumps: ```wl In[1]:= n = 2000; id = IdentityMatrix[n, SparseArray]; times = Range[0, n - 1]; temp = ConstantArray[1, n / 4]; original = Join[temp, -temp, temp, -temp] + .5 Sin[2 Pi times / n]; corrupted = original + RandomReal[{-.1, .1}, n]; ListLinePlot[{corrupted, original}] Out[1]= [image] ``` Regularize with the parameter $\lambda =2$ : ```wl In[2]:= smoothed = Fit[{id, corrupted}, "BestFitParameters", FitRegularization -> {"TotalVariation", 2.}]; In[3]:= ListLinePlot[{smoothed, original}] Out[3]= [image] ``` A smaller value of $\lambda$ gives less smoothing but the residual norm is smaller: ```wl In[4]:= lesssmoothed = Fit[{id, corrupted}, "BestFitParameters", FitRegularization -> {"TotalVariation", .1}]; In[5]:= ListLinePlot[{lesssmoothed, original}] Out[5]= [image] In[6]:= {Norm[original - smoothed], Norm[original - lesssmoothed]} Out[6]= {0.768791, 1.19395} ``` --- Use LASSO (L1) regularization to find a sparse fit (basis pursuit): ```wl In[1]:= σ = 0.05; nτ = 500; times = Subdivide[0., 1., nτ]; nω = 30; ``` Here is a signal: ```wl In[2]:= signal = (1 + .5 Sin[11 times])Sin[30 Sin[5 times]]; sdata = Transpose[{times, signal}]; ListLinePlot[sdata] Out[2]= [image] ``` The goal is to approximate the signal with just of few of the thousands of Gabor basis functions: ```wl In[3]:= Length[basis = Flatten[ Table[Exp[-(((t - τ)/σ))^2] {1, Table[{Cos[ω t], Sin[ω t]}, {ω, Drop[Subdivide[0., 150., nω], 1]}]}, {τ, Subdivide[0., 1., nτ]}]]] Out[3]= 30561 ``` With $\lambda =1.$, a fit is found using only 41 of the basis elements: ```wl In[4]:= AbsoluteTiming[sparseFit = Fit[sdata, basis, t, "BestFitParameters", FitRegularization -> {"LASSO", 1.}]] Out[4]= {3.51179, SparseArray[Automatic, {30561}, 0., {1, {{0, 45}, {{367}, {671}, {852}, {913}, {2377}, {3222}, {4195}, {4256}, {4798}, {5346}, {5890}, {6677}, {7159}, {7341}, {9456}, {9517}, {11916}, {12039}, {12411}, {13393}, {13454}, {14008}, {14375}, {15114}, {16032}, {17444}, {17505}, {18666}, {19581}, {19642}, {21043}, {21104}, {21834}, {22437}, {23349}, {23891}, {24561}, {25105}, {25770}, {26435}, {26800}, {26861}, {28976}, {29098}, {29688}}}, CompressedData["«550»"]}]} ``` The error is quite small: ```wl In[5]:= fitt = basis.sparseFit /. t -> times; ListLinePlot[Transpose[{times, signal - fitt}]] Out[5]= [image] ``` Once the important elements of the basis have been found, error can be reduced by finding the least squares fit to these elements: ```wl In[6]:= bpos = Flatten[sparseFit["NonzeroPositions"]]; sbasis = basis[[bpos]]; sbfit = Fit[sdata, sbasis, t]; ListLinePlot[Transpose[{times, signal - (sbfit /. t -> times)}], PlotRange -> All] Out[6]= [image] ``` With a smaller value $\lambda =0.1$, a fit is found using more of the basis elements: ```wl In[7]:= AbsoluteTiming[sparseFit = Fit[sdata, basis, t, "BestFitParameters", FitRegularization -> {"LASSO", 0.1}]] Out[7]= {6.92857, SparseArray[Automatic, {30561}, 0., {1, {{0, 54}, {{61}, {791}, {1282}, {1329}, {1699}, {2316}, {3100}, {4317}, {4742}, {5285}, {5957}, {6799}, {7037}, {7463}, {8112}, {8833}, {9336}, {9639}, {11855}, {11916}, {12161}, {12167} ... {17993}, {18239}, {19154}, {19886}, {20738}, {21897}, {22071}, {22444}, {23227}, {23586}, {23715}, {24622}, {25105}, {25770}, {26007}, {26679}, {26800}, {27876}, {28854}, {29871}, {29952}, {29986}, {30180}}}, CompressedData["«652»"]}]} ``` The error is even smaller: ```wl In[8]:= fitt = basis.sparseFit /. t -> times; ListLinePlot[Transpose[{times, signal - fitt}]] Out[8]= [image] ``` ### Properties & Relations (5) ``Fit`` gives the best fit function: ```wl In[1]:= Fit[Range[10], {1, x ^ 2}, x] Out[1]= 2.17582  + 0.0863422 x^2 ``` ``LinearModelFit`` allows for extraction of additional information about the fitting: ```wl In[2]:= lm = LinearModelFit[Range[10], x ^ 2, x] Out[2]= FittedModel[2.17582  + 0.0863422 x^2] ``` Extract the fitted function: ```wl In[3]:= lm["BestFit"] Out[3]= 2.17582  + 0.0863422 x^2 ``` Extract additional results and diagnostics: ```wl In[4]:= lm[{"RSquared", "FitResiduals", "ANOVATable"}] Out[4]= {0.949765, {-1.26217, -0.521193, 0.0470958, 0.4427, 0.66562, 0.715856, 0.593407, 0.298273, -0.169545, -0.810047}, | "" | "DF" | "SS" | "MS" | "F‐Statistic" | "P‐Value" | | :------ | :--- | :------ | :------- | :----------- ... | x^2 | 1 | 78.3556 | 78.3556 | 151.25 | 1.7775387117245832`*^-6 | | "Error" | 8 | 4.14443 | 0.518053 | | | | "Total" | 9 | 82.5 | | | |} ``` --- Here is some data: ```wl In[1]:= data = {{0, 1}, {1, 0}, {3, 2}, {5, 4}}; ``` This is the sum of squares error for a line `` a + b x`` : ```wl In[2]:= ss[a_, b_] = Block[{x = data[[All, 1]], y = data[[All, 2]], res}, res = y - (a + b x);res.res] Out[2]= (1 - a)^2 + (4 - a - 5 b)^2 + (2 - a - 3 b)^2 + (-a - b)^2 ``` Find the minimum symbolically: ```wl In[3]:= Minimize[ss[a, b], {a, b}] Out[3]= {(96/59), {a -> (11/59), b -> (41/59)}} In[4]:= N[a + b x /. %[[2]]] Out[4]= 0.186441  + 0.694915 x ``` These are the coefficients given by ``Fit`` : ```wl In[5]:= Fit[data, {1, x}, x] Out[5]= 0.186441  + 0.694915 x ``` The exact coefficients may be found by using ``WorkingPrecision -> Infinity`` : ```wl In[6]:= Fit[data, {1, x}, x, "BestFitParameters", WorkingPrecision -> Infinity] Out[6]= {(11/59), (41/59)} ``` This is the sum of squares error for a quadratic ``a + b x + c x ^ 2`` : ```wl In[7]:= ss[a_, b_, c_] = Block[{x = data[[All, 1]], y = data[[All, 2]], res}, res = y - (a + b x + c x ^ 2);res.res] Out[7]= (1 - a)^2 + (4 - a - 5 b - 25 c)^2 + (2 - a - 3 b - 9 c)^2 + (-a - b - c)^2 ``` Find the minimum symbolically: ```wl In[8]:= Minimize[ss[a, b, c], {a, b, c}] Out[8]= {(128/199), {a -> (135/199), b -> -(53/199), c -> (38/199)}} In[9]:= N[a + b x + c x ^ 2 /. %[[2]]] Out[9]= 0.678392  - 0.266332 x + 0.190955 x^2 ``` These are the coefficients given by ``Fit`` : ```wl In[10]:= Fit[data, {1, x, x ^ 2}, x, "BestFitParameters", WorkingPrecision -> Infinity] Out[10]= {(135/199), -(53/199), (38/199)} ``` --- When a polynomial fit is done to a high enough degree, ``Fit`` returns the interpolating polynomial: ```wl In[1]:= data = {{0, 1}, {1, 0}, {3, 2}, {5, 4}}; Fit[data, {1, x, x ^ 2, x ^ 3}, x] Out[1]= 1.  - 2.06667 x + 1.2 x^2 - 0.133333 x^3 ``` The result is consistent with that given by ``InterpolatingPolynomial`` : ```wl In[2]:= InterpolatingPolynomial[data, x] Out[2]= 1 + (-1 + ((2/3) - (2/15) (-3 + x)) (-1 + x)) x In[3]:= N[Expand[%]] Out[3]= 1.  - 2.06667 x + 1.2 x^2 - 0.133333 x^3 ``` --- ``Fit`` will use the time stamps of a ``TimeSeries`` as variables: ```wl In[1]:= ts1 = TemporalData[TimeSeries, {{{1.102448341846048, 6.331833897009389, 7.5843632999043455, 44.59969789743589, 98.76273258351092, 150.00729051003992, 644.7367359379259, 1724.3968817536402, 5903.697723461762, 19468.09340765025}}, {{0, 9, 1}}, 1, {"Continuous", 1}, {"Discrete", 1}, 1, {ResamplingMethod -> {"Interpolation", InterpolationOrder -> 1}}}, False, 10.1]; In[2]:= ts1["Times"] Out[2]= {0, 1, 2, 3, 4, 5, 6, 7, 8, 9} In[3]:= base = {1, Exp[x], Exp[2 x], Exp[x ^ 2]}; In[4]:= Fit[ts1, base, x] Out[4]= 4.42581  + 1.35733 E^x + 0.000208355 E^2 x - 3.4629357536950484`*^-32 E^x^2 ``` Rescale the time stamps and fit again: ```wl In[5]:= ts2 = TimeSeriesRescale[ts1, {.1, 2}] Out[5]= TemporalData[TimeSeries, {{{1.102448341846048, 6.331833897009389, 7.5843632999043455, 44.59969789743589, 98.76273258351092, 150.00729051003992, 644.7367359379259, 1724.3968817536402, 5903.697723461762, 19468.09340765025}}, {{0.1, 2., 0.2111111111111111}}, 1, {"Continuous", 1}, {"Discrete", 1}, 1, {ResamplingMethod -> {"Interpolation", InterpolationOrder -> 1}}}, False, 14.3] In[6]:= ts2["Times"] Out[6]= {0.1, 0.311111, 0.522222, 0.733333, 0.944444, 1.15556, 1.36667, 1.57778, 1.78889, 2.} In[7]:= Fit[ts2, base, x] Out[7]= -1385.26 + 1058.5 E^x - 428.486 E^2 x + 666.84 E^x^2 ``` Find fit for the values: ```wl In[8]:= Fit[ts1["Values"], base, x] Out[8]= 4.42581  + 0.499335 E^x + 0.0000281978 E^2 x - 1.9402117413146752`*^-40 E^x^2 ``` --- ``Fit`` acts pathwise on a multipath ``TemporalData`` : ```wl In[1]:= Fit[TemporalData[Automatic, {{{1.102448341846048, 6.331833897009389, 7.5843632999043455, 44.59969789743589, 98.76273258351092, 150.00729051003992, 644.7367359379259, 1724.3968817536402, 5903.697723461762, 19468.09340765025}, {1.1024483418 ... 7359379259, 1724.3968817536402, 5903.697723461762, 19468.09340765025}}, {{0, 9, 1}, {0.1, 2., 0.2111111111111111}}, 2, {"Continuous", 2}, {"Discrete", 2}, 1, {ResamplingMethod -> {"Interpolation", InterpolationOrder -> 1}}}, False, 10.1], {1, Exp[x], Exp[2 x]}, x] Out[1]= {-35.3494 + 1.69253 E^x + 0.0000883288 E^2 x, 5457.74  - 4804.96 E^x + 879.85 E^2 x} ``` ### Possible Issues (2) Here is some data from a random perturbation of a Gaussian: ```wl In[1]:= data = Block[{r = RandomReal[100, 100]}, Transpose[{r, Exp[-(r - 50) ^ 2 / 10] + RandomReal[.1, 100]}]]; lp = ListPlot[data, PlotRange -> All] Out[1]= [image] ``` This is a function that gives the standard basis for polynomials: ```wl In[2]:= nbasis[n_Integer][x_] := Table[x ^ k, {k, 0, n}] In[3]:= nbasis[10][x] Out[3]= {1, x, x^2, x^3, x^4, x^5, x^6, x^7, x^8, x^9, x^10} ``` Show the fits computed for successively higher-degree polynomials: ```wl In[4]:= Table[ fit = Fit[data, nbasis[n][x], x]; Show[Plot[fit, {x, 0, 100}, PlotRange -> All], lp, PlotRange -> {-.25, 1}, PlotLabel -> "n = " <> ToString[n]], {n, 10, 40, 10}] Out[4]= [image] ``` The problem is that the coefficients get very small for higher powers: ```wl In[5]:= Fit[data, nbasis[10][x], x] Out[5]= -0.203451 + 0.259948 x - 0.0740073 x^2 + 0.00935488 x^3 - 0.000626958 x^4 + 0.0000244422 x^5 - 5.808416252166985`*^-7 x^6 + 8.523541061057498`*^-9 x^7 - 7.537257251778492`*^-11 x^8 + 3.6811638504550777`*^-13 x^9 - 7.627243804704008`*^-16 x^10 ``` Giving the basis in terms of scaled and shifted values helps with this problem: ```wl In[6]:= Table[ fit = Fit[data, nbasis[n][(x - 50) / 50], x]; Show[Plot[fit, {x, 0, 100}, PlotRange -> All], lp, PlotRange -> {-.25, 1}, PlotLabel -> "n = " <> ToString[n]], {n, 10, 40, 10}] Out[6]= [image] ``` --- Reconstruct a sparse signal from compressed sensing data: ```wl In[1]:= k = 7; n = 1000; signal = Normal[SparseArray[RandomInteger[{1, n}, k] -> RandomReal[{-1, 1}, k], {n}]]; ListPlot[signal, Filling -> Axis] Out[1]= [image] ``` Compressed sensing consists of multiplying the signal with an $m\times n$ matrix with $m <<\text{\textit{$n$}}$ so only data of length $m$ has to be transmitted. If $m$ is large enough, a random matrix with independently identically distributed normal entries will satisfy the restricted isometry property, and the original signal can be recovered with very high probability: ```wl In[2]:= m = Ceiling[6 k Log[n / k]] Out[2]= 209 In[3]:= a = RandomVariate[NormalDistribution[0, 1 / m], {m, n}]; ``` The samples to send are constructed by multiplication of the signal with the matrix $a$ : ```wl In[4]:= samples = a.signal; ``` Reconstruction can be done by minimizing $\| x\| _1$ for all possible solutions of $a.x=b$ : ```wl In[5]:= AbsoluteTiming[reconstructed = NArgMin[{Norm[x, 1], a.x == samples}, x];] Out[5]= {11.3756, Null} In[6]:= Max[Abs[reconstructed - signal]] Out[6]= 1.0853510278019396`*^-8 ``` Even though the minimization is solved by linear optimization, it is relatively slow because all the constraints are equality constraints. The solution can be found much faster using basis pursuit (L1 regularization): ```wl In[7]:= AbsoluteTiming[bpsol = Fit[{a, samples}, FitRegularization -> {"L1", 0.0001}]] Out[7]= {0.0107339, SparseArray[Automatic, {1000}, 0., {1, {{0, 7}, {{35}, {189}, {348}, {681}, {808}, {902}, {922}}}, {0.47842722077532424, -0.41380898394841914, -0.0922556414811242, -0.6989351028231252, 0.9004387457771771, -0.04774504086265376, 0.16593240357071015}}]} ``` This gives the basis elements. To find the best solution, solve the linear equations corresponding to those components: ```wl In[8]:= basisPositions = bpsol["NonzeroPositions"]; values = LinearSolve[a[[All, Flatten[basisPositions]]], samples]; bpReconstructed = SparseArray[basisPositions -> values, {n}] Out[8]= SparseArray[Automatic, {1000}, 0, {1, {{0, 7}, {{35}, {189}, {348}, {681}, {808}, {902}, {922}}}, {0.4876796403830947, -0.4268590489824304, -0.10414929524049824, -0.7101969122816207, 0.9097758401442839, -0.062167722181538355, 0.17564599811089585}}] In[9]:= Max[Abs[bpReconstructed - signal]] Out[9]= 6.661338147750939`*^-16 ``` ## See Also * [`FindFit`](https://reference.wolfram.com/language/ref/FindFit.en.md) * [`LinearModelFit`](https://reference.wolfram.com/language/ref/LinearModelFit.en.md) * [`NonlinearModelFit`](https://reference.wolfram.com/language/ref/NonlinearModelFit.en.md) * [`DesignMatrix`](https://reference.wolfram.com/language/ref/DesignMatrix.en.md) * [`LeastSquares`](https://reference.wolfram.com/language/ref/LeastSquares.en.md) * [`PseudoInverse`](https://reference.wolfram.com/language/ref/PseudoInverse.en.md) * [`Interpolation`](https://reference.wolfram.com/language/ref/Interpolation.en.md) * [`InterpolatingPolynomial`](https://reference.wolfram.com/language/ref/InterpolatingPolynomial.en.md) * [`Solve`](https://reference.wolfram.com/language/ref/Solve.en.md) * [`ListPlot`](https://reference.wolfram.com/language/ref/ListPlot.en.md) * [`ConvexOptimization`](https://reference.wolfram.com/language/ref/ConvexOptimization.en.md) ## Tech Notes * [Manipulating Numerical Data](https://reference.wolfram.com/language/tutorial/NumericalOperationsOnData.en.md#25465) * [Curve Fitting](https://reference.wolfram.com/language/tutorial/NumericalOperationsOnData.en.md#19037) * [Implementation notes: Numerical and Related Functions](https://reference.wolfram.com/language/tutorial/SomeNotesOnInternalImplementation.en.md#20880) ## Related Guides * [Curve Fitting & Approximate Functions](https://reference.wolfram.com/language/guide/CurveFittingAndApproximateFunctions.en.md) * [Data Transforms and Smoothing](https://reference.wolfram.com/language/guide/DataTransformsAndSmoothing.en.md) * [Matrix-Based Minimization](https://reference.wolfram.com/language/guide/MatrixBasedMinimization.en.md) * [Numerical Data](https://reference.wolfram.com/language/guide/NumericalData.en.md) ## Related Links * [NKS\|Online](http://www.wolframscience.com/nks/search/?q=Fit) * [A New Kind of Science](http://www.wolframscience.com/nks/) ## History * Introduced in 1988 (1.0) \| [Updated in 2019 (12.0)](https://reference.wolfram.com/language/guide/SummaryOfNewFeaturesIn120.en.md)