5
$\begingroup$

I'm trying to fit the following data using Mathematica, but unfortunately I'm not getting any decent results.

data = {{1026689495,0.04346642},{1032948957,0.07106946},{1033037820,0.07173154},{1033175985,0.07276095},{1033544110,0.07550369},{1033794716,0.07737085},{1035382618,0.08920164},{1035473824,0.08988118},{1036246565,0.09563855},{1036526395,0.09772345},{1046404258,0.21763852},{1050599907,0.31031713},{1053371976,0.39293122},{1058188482,0.57294858},{1062404620,0.82424998},{1064270054,0.93543851},{1002403160,0.00592598},{978672841,0.00081391}};

I'm assuming the following model:

model = a Exp[k t]

I tried using FindFit:

fit = FindFit[data, model, {a, k}, t]

Which errors with

Failed to converge to the requested accuracy or precision within 100 iterations.

I tried NonlinearModelFit too but I can't get it to work.

In the end I tried using Wolfram Alpha which gives decent results but I can't fit all of my data into the textfield and it doesn't work with WolframAlpha Pro.

4.15092×10^-39 e^(0.00830101 x)

How would one solve this using Mathematica? Wolfram Alpha Fit

$\endgroup$
6
  • 1
    $\begingroup$ Did you try changing the MaxIterations to something larger than 100 ? $\endgroup$
    – Sektor
    Jul 6, 2015 at 15:09
  • $\begingroup$ Yes i even tried increasing it to 100 000 and the Results are still far far worse than WolframAlpha. :/ $\endgroup$ Jul 6, 2015 at 15:10
  • 2
    $\begingroup$ Well, with such extreme values you should either shift the plot points, or specify initial values. Do you know what a and k should roughly be? $\endgroup$
    – Feyre
    Jul 6, 2015 at 15:15
  • $\begingroup$ Thank you. I'll try shifting the plot points. I already tried specifying initial values(the ones from WA), but the results were still not good. Unfortunately i dont know what a and k should roughly be(again - except the WA guesses) $\endgroup$ Jul 6, 2015 at 15:21
  • 1
    $\begingroup$ Do you have all you need now, because it's hard to guess without knowing where the data physically comes from? $\endgroup$
    – Feyre
    Jul 6, 2015 at 15:42

3 Answers 3

11
$\begingroup$

If you're pretty sure it's an exponential, you can always take the logarithm of the data and do a linear fit to that:

logdata = {#[[1]], Log[#[[2]]]} & /@ data;
FindFit[logdata, c x + d, {c, d}, x]

(* {c -> 8.2386*10^-8, d -> -87.7291} *)

Note that c and d are related to the original parameters by $c = k$ and $d = \ln a$. This means that the correct value for $a$ in your original fit was about $e^{-87} \approx 8 \times 10^{-39}$, which explains why Mathematica was having trouble finding a numerical answer for it.

For reference, here's what your data looks like on a logarithmic plot, with the fit superposed:

Show[ListLogPlot[data], Plot[c x + d /. fit, {x, 9.8*10^8, 1.07*10^9}, PlotStyle -> {Red}]]

enter image description here

And here are the "fractional residuals", defined as the ratio between the best fit and the actual data points:

{t, y} = Transpose[data]
ListLogPlot[Transpose[{t, y/((Exp[c # + d] /. fit) & /@ t)}], AxesOrigin -> {Automatic, 1}, Filling -> Axis]

enter image description here

Note that if we're fitting logarithmically scaled data with a linear fit, then we're effectively trying to minimize the sum of the squares of the ratios between the data points and the fit values. This would be the best fit if we expect the percent errors of the data points to be normally distributed (and all the same), rather than the absolute errors. If the data have different error bars than this, then you might want to pursue a different technique.

$\endgroup$
10
  • 1
    $\begingroup$ Where are you getting the number that a "should be"? $\endgroup$ Jul 6, 2015 at 16:34
  • 3
    $\begingroup$ Quite possibly, yes. From my edit, you can see that your data vary from about 1-2% from a pure exponential. Mathematica finds the best fit assuming that all the data points are equally uncertain. If you happen to know that the first data point is more accurate than all the rest, then there are ways to make the fit get closer to that point, but at the cost of moving the fit away from the other 17 points (on average.) $\endgroup$ Jul 6, 2015 at 17:28
  • 3
    $\begingroup$ @Michael Seifert. I'd go farther and argue it's completely unreasonable to obtain one more digit of precision. In fact, because of what you show with the residuals (observations typically within about 1 to 2% of the predicted value) one should consider at most 3 to 4 digits to the right of the decimal to be accurate (or find a better model). $\endgroup$
    – JimB
    Jul 6, 2015 at 19:01
  • 2
    $\begingroup$ @MichaelSeifert I would argue here that doing a nonlinear least squares fit on the original data is not the same as linear least squares fit in log-log space. If you assume that the original data satisfies the assumptions of least squares then taking the log of the data will violate those assumptions and LS will give wrong/biased results. One could take the linear LS estimates and use those as initial conditions for the correct nonlinear LS fit. $\endgroup$ Jul 7, 2015 at 0:27
  • 2
    $\begingroup$ @MichaelSeifert I have held the same impression myself and god only knows how many such analyses I have done. But if you ponder for a moment you will see that weighted NLLS and unweighted LLS cannot be equal. In the weighted case, you are fitting (wx,wy) where w includes weights. In the nonweighted case, you are fitting (log(x),log(y)). The difference is that w is a constant and it is the SAME on both x and y. Taking the log is equivalent to different "weights" on both sides unless x and y are identical (in which case the fitting is unnecessary because y=x will fit the data). $\endgroup$ Jul 7, 2015 at 3:42
6
$\begingroup$

The original issue of a lack of convergence is because of the very large values of the independent variable and the use of the default starting values (and my understanding is that the default starting value for all parameters is 1.0).

This can be fixed by standardizing the independent variable. (This is not a bad practice for just about any regression analysis - but you have to remember to make the conversion back to the original units.)

(* Data *)
data = {{1026689495, 0.04346642}, {1032948957, 
    0.07106946}, {1033037820, 0.07173154}, {1033175985, 
    0.07276095}, {1033544110, 0.07550369}, {1033794716, 
    0.07737085}, {1035382618, 0.08920164}, {1035473824, 
    0.08988118}, {1036246565, 0.09563855}, {1036526395, 
    0.09772345}, {1046404258, 0.21763852}, {1050599907, 
    0.31031713}, {1053371976, 0.39293122}, {1058188482, 
    0.57294858}, {1062404620, 0.82424998}, {1064270054, 
    0.93543851}, {1002403160, 0.00592598}, {978672841, 0.00081391}};

(* Model form being fit *)
model = a Exp[k t];

Using default starting values with the original data...

  FindFit[data, model, {a, k}, t]
  FindFit[data, model, {{a, 1}, {k, 1}}, t]

results in

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

Clearly not a good enough answer. Choosing better starting values helps but sometimes only if the number of iterations is increased from the default.

FindFit[data, model, {{a, 10^(-38)}, {k, 10^(-7)}}, t]
FindFit[data, model, {{a, 10^(-38)}, {k, 10^(-7)}}, t, 
 MaxIterations -> 200]

resulting in

FindFit::cvmit: Failed to converge to the requested accuracy or precision within 100 iterations. >>
{a -> 1.28015*10^-38, k -> 8.19291*10^-8}
{a -> 1.72511*10^-38, k -> 8.16492*10^-8}

Had we standardized the independent variable there would have been no need to increase the iterations or have a good guess as to the starting values:

data2 = data;
data2[[All, 1]] = Standardize[data[[All, 1]]];
sol = FindFit[data2, model, {a, k}, t];
xbar = N[Mean[data[[All, 1]]]];
sigma = N[StandardDeviation[data[[All, 1]]]];
{a Exp[-k xbar/sigma], k/sigma} /. sol

with output

{1.72511*10^-38, 8.16492*10^-8}

Having said all of the above for this particular dataset a better approach that more closely matches the residual error structure is the taking of the logs of the dependent variable and performing a linear regression as Michael Seifert did (but I'd use LinearModelFit rather than NonlinearModelFit - however, in this case the results are equivalent). Using LinearModelFit (or NonlinearModelFit) provides a whole lot more information about the fit than FindFit.

In doing so one would see that requesting more decimal places in the predictions does not provide a better fit and that only 3 or maybe 4 digits to the right of the decimal are warranted given the quality of the fit. And with only 18 data points a more complex model with maybe a better fit to the observed data is not justifiable.

$\endgroup$
1
  • $\begingroup$ +1 Useful for other things, I think. $\endgroup$
    – Carl
    Aug 6, 2022 at 17:56
5
$\begingroup$

First things first, let's shift the data to the origin.

data = data - ConstantArray[{10*^8, 0}, 18];

Try some nice initial values:

nlm = NonlinearModelFit[data, a Exp[k t], {{a, 0.01}, {k, 0.0000001}},  t];
nlm["BestFitParameters"]
Show[ListPlot[data], Plot[nlm[t], {t, 0, 7*^7}], Frame -> True]
(* {a -> 0.00497293, k -> 8.16492*10^-8} *)

Mathematica graphics

And that should do it.

$\endgroup$
2
  • $\begingroup$ Thank you, this gives close values. The actual values are still not precise enough(should be: 0.0434664 nlm: 0.0439558), is there any way to increase the precision? Increasing MaxIterations did not help. $\endgroup$ Jul 6, 2015 at 15:50
  • 2
    $\begingroup$ Well, the origin could be off. $\endgroup$
    – Feyre
    Jul 6, 2015 at 16:02

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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