How to improve working precision?

Question

I'm trying to fit the array data as follow:

NonlinearModelFit[
  data, 
  s*t^r + (A/(w*(Pi/2)^0.5))*Exp[-2*((t - u)/w)^2], 
  {s, r, A, w, u}, 
  t];

But when I run it, it gives the following problem:

Exp[] is too small to represent as a normalized machine number; precision may be lost.

Do you know how to solve it?

Thank you!

Here the code:

data={{9331200, 2.25014*10^-14}, {9590400, 2.94221*10^-14}, {12441600, 
  2.25158*10^-14}, {13219200, 1.88385*10^-14}, {13564800, 
  8.74689*10^-15}, {13824000, 1.65873*10^-14}, {13910400, 
  1.91969*10^-14}, {14083200, 1.85371*10^-14}, {14256000, 
  1.21808*10^-14}, {22464000, 1.55479*10^-14}};

errors={2.35061*10^-15, 4.03579*10^-15, 2.57128*10^-15, 3.17395*10^-15, 
 4.44869*10^-15, 3.86327*10^-15, 3.73115*10^-15, 1.55858*10^-15, 
 4.62527*10^-15, 1.50732*10^-15};

dataPlot = ListLogLogPlot[data, PlotRange -> Full];

(*First Fit with Powerlaw, it works*)
Pow = NonlinearModelFit[data, k*t^h, {k,h}, t, 
   MaxIterations -> 1000, Weights -> 1/errors^2];
PowPlot = 
  LogLogPlot[Pow[t], {t, 100*24*3600, 260*24*3600}, PlotStyle -> Red, 
   PlotRange -> Full];
"Plot powerlaw:"
Show[{PowPlot, dataPlot}, PlotRange -> Full]

Pow["ParameterTable"]

chi1r = Pow["EstimatedVariance"]


(*Second Fit with Powerlaw+Gaussian, where it gives that problem*)
Pgauss = NonlinearModelFit[data, 
   s*t^r + (A/(w*(Pi/2)^(1/2)))*Exp[-2*((t - u)/w)^2], {s, r, A, w, 
    u}, t, Weights -> 1/errors^2, MaxIterations -> 1000, 
   WorkingPrecision -> 50, Method -> NMinimize];

PgaussPlot = 
  LogLogPlot[Pgauss[t], {t, 100*24*3600, 260*24*3600}, PlotStyle -> Red];

"Plot powerlaw+gaussiana:"
Show[{PgaussPlot, dataPlot}]

Pgauss["ParameterTable"]
chi2r = Pgauss["EstimatedVariance"]

Thank you very much!

Answer 1

Higher precision is not needed. One just needs to scale the data so that numerical instabilities don't occur. This is especially the case here where the scales of the predictor and response variables are so far apart.

For both models multiplying the predictor and response by appropriate constants don't change the underlying model. (Sometimes subtracting a constant is what you want which doesn't change the underlying model but that's not the case here.)

(* Original data *)
data = {{9331200, 2.25014*10^-14}, {9590400, 2.94221*10^-14}, {12441600, 2.25158*10^-14},
   {13219200, 1.88385*10^-14}, {13564800, 8.74689*10^-15}, {13824000,  1.65873*10^-14},
   {13910400, 1.91969*10^-14}, {14083200, 1.85371*10^-14}, {14256000,  1.21808*10^-14},
   {22464000, 1.55479*10^-14}};
errors = {2.35061*10^-15, 4.03579*10^-15, 2.57128*10^-15, 3.17395*10^-15, 4.44869*10^-15,
   3.86327*10^-15, 3.73115*10^-15, 1.55858*10^-15, 4.62527*10^-15, 1.50732*10^-15};

(* Scale data *)
errors2 = errors 10^15;
data2 = data;
data2[[All, 1]] = data[[All, 1]]/10^7;
data2[[All, 2]] = data[[All, 2]] 10^14;

(* Power law fit *)
Pow = NonlinearModelFit[data2, k*t^h, {k, h}, t, Weights -> 1/errors2^2];
Pow["ParameterTable"]

(* More complicated model *)
Pgauss = NonlinearModelFit[data2, s*t^r + (A/(w (Pi/2)^(1/2)))*Exp[-2*((t - u)/w)^2],
   {s, r, A, w, u}, t, Weights -> 1/errors2^2];

(* Results *)
Pgauss["ParameterTable"]

 Show[ListPlot[data2],
   Plot[{Pow[t], Pgauss[t]}, {t, Min[data2[[All, 1]]], Max[data2[[All, 1]]]}, PlotRange -> All,
   PlotLegends -> {"Pow model", "Pgauss model"}], PlotRange -> All]

The Pgauss fit looks pretty implausible as an explanation. (The fit for the Pgauss model is fine. It's the consideration of using it to describe the relationship is what is implausible.) The $AIC_c$ statistics bear that out:

Pow["AICc"]
(* 15.9845 *)
Pgauss["AICc"]
(* 41.1028 *)

as the $AIC_c$ statistic for the Pow model is far lower than that of the Pgauss model.