4
$\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$
  • 1
    $\begingroup$ Did you try changing the MaxIterations to something larger than 100 ? $\endgroup$ – Sektor Jul 6 '15 at 15:09
  • $\begingroup$ Yes i even tried increasing it to 100 000 and the Results are still far far worse than WolframAlpha. :/ $\endgroup$ – Marcel Kiesel Jul 6 '15 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 '15 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$ – Marcel Kiesel Jul 6 '15 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 '15 at 15:42
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$
  • 1
    $\begingroup$ Where are you getting the number that a "should be"? $\endgroup$ – Michael Seifert Jul 6 '15 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$ – Michael Seifert Jul 6 '15 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 '15 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$ – Fixed Point Jul 7 '15 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$ – Fixed Point Jul 7 '15 at 3:42
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$
  • $\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$ – Marcel Kiesel Jul 6 '15 at 15:50
  • 2
    $\begingroup$ Well, the origin could be off. $\endgroup$ – Feyre Jul 6 '15 at 16:02
5
$\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$

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.