35
$\begingroup$

Mathematica has numerous functions designed to, or capable of, fitting known functions, and finding unknown functions to match data sets. What are some common issues that come with finding those fits?

$\endgroup$
  • 2
    $\begingroup$ You may want to mention EstimatedDistribution and FindDistributionParameters along with a few words about fitting curves vs fitting distributions. $\endgroup$ – Szabolcs Mar 2 '17 at 16:57
40
$\begingroup$

First, let's enumerate some of the functions:

  • Fit[] is the simplest of fitting functions. It has been introduced in v5, and hasn't been updated since v6 (as of v11). It finds a least-squares fit as a linear combination of functions. As such it can only be used for simple functions. Fit[] can fit $c\sin{x}$ but not $\sin{c x}$, where c is an unknown parameter for Fit[] to find.

  • FindFit[] has been introduced in v5, and hasn't been updated since v7 (as of v11). It provides a simple and clear way of finding parameters of a known function to a dataset. FindFit[] also finds a simple least-squares fit, but can be used for more complicated functions, including Trigonometric functions.

  • LinearModelFit[] and NonlinearModelFit[] have been introduced in v7 and have been continuously updated since, which are related to Fit[] and FindFit[] respectively.

    These functions construct a fit with a lot of additional information which can be extracted.

    In the same vain as LinearModelFit[]there are some additional symbols such as LogitModelFit and ProbitModelFit which fit data to predefined models (as opposed to a linear model).

  • FindGeneratingFunction[] is a symbol that attempts to generate a function to a dataset. The symbol will attempt to find a function that matches the data points exactly.

  • FindFormula[] is a symbol which attempts to find a formula to non-exact data. FindFormula[] is currently unreliable.

  • FindSequenceFunction[] is a symbol which tries to exactly match known series to a dataset.

  • Interpolation[] generates a piecewise function which interpolates between the points it is fed. It doesn't return an analytical function.

  • NMinimize[] minimizes a function with respect to a variable/ variables. This function can be used to create a fitting function.

Fitting: Poor or no results

  • 0. Insufficient or erroneous data points

    Though this should be obvious, it is nonetheless a common problem. As a general rule, in order to find n parameters, you need n data points. However, considering many of these functions are not exact, usually many more data points are required. In general, more data points is better. Additionally data points should be fed in the form {{x1, y1}, {x2, y2}, {x2, y3}, ..., {xn, yn}}. If no the data is in the form {y1, y2, y3, ..., yn} then the x-values will be taken to be Range[n].


  1. Initial Values

    Often times the estimates that Mathematica produces are poor, especially when the real values are far from 1 or there are large order differences between the variables.

    Consider:

    data = {{7, 17.93}, {14, 36.36}, {21, 67.76}, {28, 98.1}, {35, 
        131}, {42, 169.5}, {49, 205.5}, {56, 228.3}, {63, 247.1}, {70, 
        250.5}, {77, 253.8}, {84, 254.5}};
    model = c/(1 + ((c - n)/n)*E^(-r*t));
    fit = FindFit[data, model, {c, n, r}, t];
    

    The result is quite poor:

    Show[Plot[model /. fit, {t, 7, 84}], ListPlot@data]
    

    enter image description here

    Initial estimates of the function can be used. The list of variables has to be augmented with the estimates: {{c, 256}, {n, 9}, r} instead of {c, n, r}. The result is significantly better:

    fit = FindFit[data, model, {{c, 256}, {n, 9}, r}, t]
    

    enter image description here

  • 1.1

    Supposing difficulties with finding the initial values, NonlinearModelFit[] can be used to programmatically find good values.

    data = Table[{i, 84.95 Sin[2 i] + 50}, {i, 1, 20}];
    model = {a Sin[b x + c] + d, -π <= c <= π};
    nlm = Quiet@Table[
       NonlinearModelFit[
        data, model, {a, {b, i}, c, 
         d}, x, MaxIterations -> 1000], {i, 0, 4, 0.05}];
    

    This will generate a list of fits, which can be tested.

    ListPlot[Norm[nlm[[#]]["FitResiduals"]] & /@ Range@81, Joined -> True]
    

    enter image description here

    As you can see, there are sharp spikes where the fit converges, finding the first and acquiring the fit:

    fit = Normal@nlm[[Position[
     min = Norm[nlm[[#]]["FitResiduals"]] & /@ Range@81, Min[min]][[1,
      1]]]]
    

    50. + 84.95 Sin[3.52466*10^-17 + 2. x]

    Show[Plot[fit, {x, 1, 20}], ListPlot@data]
    

    enter image description here

  • 1.2 Overflow in fit

    Exponential data, and other data that has big differences in y-values of data can create overflow in the fit.

    data = {{0, 6.51}, {100, 0.77}, {200, 0.306}, {400, 0.0476}, {700, 0.004}};
    model = a Exp[-b x];
    nlm = FindFit[data, model, {a, b}, x]
    

    General::ovfl: Overflow occurred in computation. >>

    Using method "NMinimize" prevents overflow.

    FindFit[data, model, {a, b}, x, Method -> "NMinimize"]
    

    {a -> 6.50674, b -> 0.0206953}

  • 1.3 Weighted data results in bad fit

    Assuming the above data.

    Using NonlinearModelFit[] you can add error data.

    One might assume that in some physical models, the error scales linearly with the magnitude of the data. Take into account that weights are rated quadratically. With limited data, using square root scaled weights is usually the right rate.

    errors = data[[All, 2]];
    nlme = NonlinearModelFit[data, model, {a, b}, x, 
      Method -> "NMinimize", Weights -> 1/errors^2]; nlme2 = 
     NonlinearModelFit[data, model, {a, b}, x, Method -> "NMinimize", 
      Weights -> 1/errors2^2];
    nlm = NonlinearModelFit[data, model, {a, b}, x, Method -> "NMinimize"];
    Show[ListPlot[data, PlotStyle -> {PointSize[0.02], Darker@Red}], 
     Plot[{nlm[x], nlme[x], nlme2[x]}, {x, 0, 700}, PlotRange -> Full], 
     PlotRange -> All]
    

    enter image description here

    Where blue is unweighted, yellow is linearly weighted, and green inverse quadratically weighted.

  • 1.4 Exponential Data

    Consider the following model:

    data = {{1032948957, 0.0710695}, {1033175985, 0.072761}, {1033794716, 
      0.0773709}, {1035473824, 0.0898812}, {1036526395, 
      0.0977235}, {1050599907, 0.310317}, {1058188482, 
      0.572949}, {1064270054, 0.935439}, {978672841, 0.00081391}};
    model = a Exp[k t];
    fit = FindFit[data, model, {a, k}, t]
    

    {a -> 0., k -> 1.}

    The result is obviously wrong. The data is presented in a difficult way, the x-values are far higher than the y-values, and the data follows an exponential curve, so initial values might not do the trick.

    One efficient way of solving this, is by taking the log of the y-data and the fit:

    logdata = {#[[1]], Log[#[[2]]]} & /@ data;
    fit = FindFit[logdata, c x + d, {c, d}, x];
    Show[Plot[c x + d /. fit, {x, 1.03*^9, 1.07*10^9}, PlotStyle -> {Red}], 
     ListLogPlot[data, PlotStyle -> {Purple, PointSize[0.02]}]]
    

    enter image description here

    Note that this is a linear regression, and as such improperly weights the potential errors in the data points.

    We can add weights based on the slope of an InterpolatingFunction constructed from the data.

    slopes = D[Interpolation[data][x], x] /. x -> data[[#, 1]] & /@ 
      Range@9;
    

    And then fit with NonlinearModelfit[]:

    fitmod = NonlinearModelFit[logdata, c x + d, {c, d}, x, 
      Weights -> 1/slopes^2]
    

    The (log-)difference in the fits is small, but can be relevant.

    Plot[Normal@fitmod - (c x + d) /. fit, {x, 1.03*^9, 1.07*10^9}]
    

    enter image description here


  1. Possibly complex functions

    Consider the data set

    data = ToExpression@Import["http://pastebin.com/raw.php?i=eSXQuV62"];
    

    Trying initial values:

    model = G/(2/3*x)*(1 - (4.13/G*(Log[k/(x*r)] + G/4.13))^(2/3));
    fit = FindFit[data, fun, {{k, 50}, {x, 0.2}, {G, 80}}, r];
    

    Not only does the fitting result in a warning, but the functions cannot be plotted at all. This is because unless both x > 0 && k > 0, the model returns a complex result. As a result, functions with possible complex results need to be restrained.

    fit = Quiet@FindFit[
      data, {fun, k > 0, x > 0}, {{k, 0.5}, {x, 0.01}, {G, 80}}, r];
    Show[Plot[Re[model /. fit], {r, -10, 15000}], ListPlot@data]
    

    enter image description here

    This problem can most commonly occur with any models that have roots. Suppose some function f[x] with unknown constant c, which is real at unknown intervals

    f[x_, c_]:= ...
    

    The domain at which the function is real can be found by

    dom = Reduce[And @@ Thread[f[#, b] & /@ data[[All, 1]] > 0], {b}, Reals];

where data[[All, 1]] is the dataset's x-values.


  1. Complex Data

    • 3.1 Complex variables

    Suppose you need to fit data in the complex plane:

    data = {{5.0119`*^6, 4.3675` + 0.45799` I}, {3.9811`*^6,  4.3315` + 0.48469` I},
    {3.1623`*^6, 4.339` + 0.52722` I}, {2.5119`*^6,4.359` + 0.57605` I}, {1.9953`*^6, 4.409` + 0.63945` I},
    {1.5849`*^6,4.4774` + 0.70803` I}, {1.2589`*^6, 4.5612` + 0.78699` I}, {1.`*^6, 4.6626` + 0.87252` I}}
    model = a + b/(1 + I x c);
    

    Then

    fit = FindFit[comp, model, {a, b, c}, x];
    

    won't work, and will return the error that

    The function value ... is not a list of real numbers with dimensions

    The particular issue is the c which is paired with the I. One way of solving this is predefining c.

    model2[c_] := a + b/(1 + I x c);
    With[{c = 7*^-6}, FindFit[data, model2[c], {a, b}, x]]
    

    {a -> 4.20685 + 0.378959 I, b -> -3.17118 + 3.49223 I}

    enter image description here

    A slight modification of subsection 1.1 can be used to find an optimal value for c.

  • 3.2 Real variables

    Alternatively, a fit symbol can be constructed as follows:

    model[x_] = a + b/(1 + I x c)
    s = NMinimize[{Total[Norm[model[#[[1]]] - #[[2]]]^2 & /@ data], 
       0 < a < 10 && 0 < b < 10 && -10 < c < 0}, {a, b, c}, 
      MaxIterations -> 10000]
    

    {0.2686, {a -> 4.43805, b -> 9.44122, c -> -8.53354*10^-6}}

    enter image description here

    Again, with the limited sample data, initial domains as specified here are necessary for a good result.


  1. Implicit Functions

    Suppose you have data which has to be fitted to the form

    modelimp = a y + b Log[y] == x;
    data = {{1, 1}, {2, 1.4}, {3, 1.8}, {4, 2.4}, {5, 2.9}};
    

    The most straightforward way to solve this is to find the function in form of y =.

    model = y/.Quiet@Solve[modelimp, y][[1]];
    fit = FindFit[data, model, {a, b}, x]
    

    {a -> 0.986636, b -> 1.96879}

    Alternatively, if this isn't possible, FindRoot[] can be used to replace the model inside the fit symbol.

    fitfunc[a_?NumericQ, b_?NumericQ, x_?NumericQ] := 
      y /. FindRoot[a y + b Log[y] == x, {y, 1.}];
    fit = FindFit[data, fitfunc[a, b, x], {a, b}, x];
    Show[ListPlot[data], Plot[model /. fit, {x, 1, 5}]]
    

    enter image description here

    However, an implicit function requires a different calculation to minimize the residuals in both the x- and y- directions. To do this, you can use NArgMin to minimize the norm of the implicit function.

    {a, b} = NArgMin[
      Norm[Function[{x, y}, c[1] y + c[2] Log[y] - x] @@@ data], {c[1], 
       c[2]}]
    

    {0.990983, 1.95206}

    As you can see, this gives slightly different results for a and b.


  1. FindGeneratingFunction[] finds the wrong result or no result.

    FindGeneratingFunction[] looks for functions that match its dataspace exactly. However, if there are few datapoints, or the data points follow a logical pattern, there may be many results possible. This is especially true for infinitely repeating series.

    At other times no solution is found at all.

    Specifying , FunctionSpace -> will tell FindGeneratingFunction[] where to look. As any such series can be generated by an n-dimensional polynomial where n the length of the dataset. Polynomials of length n are often not given as a solution. Specifying , FunctionSpace -> "Polynomial" will force the simplest possible polynomial to be given as the answer, regardless of its order.


  1. Problems with prepending FindGeneratingFunction[]

    Consider

    FindGeneratingFunction[{0, 0, 0, 0, 2, 2, 6, 6, 22, 22, 86, 86, 342, 
      342, 1366, 1366, 5462, 5462}, t]
    

    (2^(1 + t) t^4 (1 - 2 t^2))/(1 - t - 4 t^2 + 4 t^3)

    The result should logically be the same as

    t^4 FindGeneratingFunction[data = {2, 2, 6, 6, 22, 22, 86, 86, 342, 342, 
       1366, 1366, 5462, 5462}, t]
    

    -((2 t^4 (-1 + 2 t^2))/(1 - t - 4 t^2 + 4 t^3))

    However, the former is incorrect. The latter is a proper workaround, programmatically this is

    t^(n = Length@Position[data, 0]) FindGeneratingFunction[Drop[data, n], t]
    

    -((2 t^4 (-1 + 2 t^2))/(1 - t - 4 t^2 + 4 t^3))

$\endgroup$
  • $\begingroup$ Hi; Nice idea to pull the shifting t^4 out and shift the data. +1. $\endgroup$ – bobbym Mar 2 '17 at 17:14
  • $\begingroup$ +1. Great resource. Wouldn‘t GeneralizedLinearModelFit be worth a note as well? $\endgroup$ – gwr Mar 2 '17 at 22:48
  • $\begingroup$ @Kuba 0 is essentially for extremely simple mistakes, but I'm not married to it. $\endgroup$ – Feyre Mar 3 '17 at 13:13
  • $\begingroup$ @Feyre let me know what do you think abut the eidt, feel free to revert it of course. $\endgroup$ – Kuba Mar 3 '17 at 13:26
11
$\begingroup$

@Feyre's answer is excellent. I would add the following two items:

1. Overparameterizing a model.

When a model has parameters that are not estimable, errors and convergence problems arise. For example a model with parameters $a$, $b$, and $c$ in the following form

$$y=a*b + c*x + error$$

does not allow one to estimate both $a$ and $b$. The product $a*b$ can be estimated but not the individual parameters. There are related issues when some of the parameter estimators are highly correlated. Looking at the "CorrelationMatrix" option for some of the fitting techniques is helpful for seeing if overparameterization is an issue.

2. Not looking at the residuals.

Most of the model fitting techniques currently assume a constant variance across all of the predictor values. Weighting is one way to address that issue but even with a well thought out weighting scheme, one needs to check on the residuals to look for lack of fit either in the prediction equation or the assumption of constant variance. Plotting the histogram of residuals and the residuals vs. the predicted values should be a matter of course to check for undesirable patterns and deviations from assumptions.

Most models have a "fixed" portion and a "random" portion with both being essential to a justifiable fit.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.