Message Boards Message Boards

0
|
7146 Views
|
5 Replies
|
4 Total Likes
View groups...
Share
Share this post:

NonlinearModelFit: how to improve accuracy, enforce constraints?

Posted 2 years ago

Dear community,

I ran into some problems with NonlinearModelFit when trying to fit experimental data. The fit does not reach a satisfying closeness and I get a few wierd error messages. After days of online search I decided to ask here.

I'll take you through my code:

At first, the data series for time and heatflow are imported (the file is attached to this post), the initial heat flow and time, and the equation to which I want to fit are defined

rawdata = SetPrecision[Import["Raw_heatflow_data.txt"], 20];
time = rawdata[[1]];
heatflow = rawdata[[2]];
ListPlot[Thread[{rawdata[[1]], rawdata[[2]]}]]
t0 = time[[1]];
q0 = heatflow[[1]];
data = Thread[{time, heatflow}];
model = (q0 - (a + b*(t - t0)))*Exp[-(t - t0)/tau] + (a + b*(t - t0));

The precision of all numbers can be verified here as 20. This will be important to an error message below.

Precision[model]
Precision[time]
Precision[heatflow]

When doing the fit, I applied any options that I found online which might potentially resolve the issues:

fit = NonlinearModelFit[
   data, {model, {10^-6`20 < a < 5*10^-3`20, -1.1*10^-9`20 < b < -1*10^-19`20, 
   10`20 < tau < 1200`20}}, {{a, 1*10^-6`20}, {b, -1*10^-9`20}, {tau, 650`20}}, t, 
   WorkingPrecision -> 20, MaxIterations -> 100, AccuracyGoal -> 15, PrecisionGoal -> 15, 
   ConfidenceLevel -> 95/100, Tolerance -> 10^-50, StepMonitor :> Print["Step to a=", a, ", b=", b, " and T=", tau]];

The resulting fit is somewhat close to the data, but when a colleague did the same fitting in Matlab, he got a much closer, almost perfect fit. Are there any options in Mathematica to improve the accuracy of the fit? I'm thinking of options like a smaller iteration step size or limiting the fit residuals to a lower value.

Also, do you have recommendations what I should do about AccuracyGoal, PrecisionGoal, ConfidenceLevel and Tolerance, as I don't completely understand what they are doing in the fitting process?

Besides these general questions, there have been some more issues with my code. As visible through the step monitor, b reaches a value of -4.96648*10^-9 at the last step which is below the lower constraint. Hence, it is not followed. Is there a way to strictly enforce the constraints at every step?

Also, two error messages pop up:

NonlinearModelFit::precw: The precision of the data and model function (MachinePrecision) is less than the specified WorkingPrecision (20).

Why is the precision of my data and model reset to MachinePrecision, and not 20 anymore, as above? How can I prevent this issue? I even wrote the `20 for all the starting values and constraints, and it still does not resolve the error...

NonlinearModelFit::sdir: Search direction has become too small.

What does this error mean and how does it affect my calculation?

Best regards

Roland

Attachments:
POSTED BY: Roland Käser
5 Replies
Posted 2 years ago

The issue seems to be that the "best" fit involves a and b set to their lower bounds and NonlinearModelFit just isn't reaching those lower values. But because it then appears obvious that the lower bounds for those two parameters are the values that maximize the likelihood, setting those two parameters to their lower bounds is necessary.

t0 = Rationalize[time[[1]], 0];
q0 = Rationalize[heatflow[[1]], 0];
data = Rationalize[Thread[{time, heatflow}], 0];
model = (q0 - (a + b (t - t0)))*Exp[-(t - t0)/tau] + (a + b*(t - t0)) /. b -> 11 10^-10 /.  a -> 10^-6;
fit = NonlinearModelFit[data, {model, {10 < tau < 1200}}, {{tau, 658}}, t];
fit["BestFitParameters"]
(* {tau -> 660.9437814041655`} *)
mseMathematica = Mean[(data[[All, 2]] - (model /. t -> data[[All, 1]] /.fit["BestFitParameters"]))^2]
(* 3.0509905090105965`*^-10 *)
mseMatlab = Mean[(data[[All, 2]] - (model /. t -> data[[All, 1]] /. {a -> 0.0000010000000222047, 
        b -> -1.09997779553949 10^-9, tau -> 657.977318176476}))^2]
(* 3.0950005923166116`*^-10 *)

One obtains slightly better fits (using the mean square error) with Mathematica than Matlab.

As Neil Signer mentioned if you want better fits you need a different model or loosen up on the parameter restrictions. Also, $R^2$ is not a great measure of goodness of fit especially if the values are near 1 as in this case. Consider using fit["EstimatedVariance"] or the square root of that value instead.

POSTED BY: Jim Baldwin
Posted 2 years ago

Looks like the raw data in the CSV file has an issue

ListPlot[rawdata[[1]],
 PlotRange -> All,
 ImageSize -> 500,
 PlotTheme -> "Detailed"]

enter image description here

POSTED BY: Rohit Namjoshi

Dear Rohit,

That "step" is on purpose as I picked two separate time segments for fitting the model.

Best regards,

Roland

POSTED BY: Roland Käser

Roland,

The problem is that your data formatting is messed up. I can't even read the data file you provided with the code you provided. By editing the data file a bit (I broke it into two files because Import was reading it as strings). I was able to read it. Next, you have all these options, constraints and precision specifications that you simply do not need. Your data does not need 20 digits of precision. For example, you have time values in 1 second increments with 20 digits of precision! Mathematica does an excellent fit out of the box.

time = Get["Raw_heatflow_data.txt"];
heatflow = Get["Raw_heatflow_data2.txt"];
ListPlot[Thread[{time, heatflow}]]
t0 = time[[1]];
q0 = heatflow[[1]];
data = Thread[{time, heatflow}];
model = (q0 - (a + b*(t - t0)))*Exp[-(t - t0)/tau] + (a + b*(t - t0));
fit = NonlinearModelFit[
   data, {model}, {{a, 1*10^-6`20}, {b, -1*10^-9`20}, {tau, 650`20}}, 
   t];
Show[Plot[fit[x], {x, time[[1]], time[[-1]]}, PlotStyle -> Red], 
 ListPlot[data], Frame -> True]

Without changing your model, you will unlikely get a better fit. How does this compare to the Matlab Result?

In[20]:= fit["BestFitParameters"]

Out[20]= {a -> -0.000075701139127212902371, 
 b -> 1.7776365345108173411*10^-8, tau -> 615.11694898053211394}

Regards,

Neil

enter image description here

POSTED BY: Neil Singer

Dear Neil,

Thank you very much for your response. Sorry about the txt file. I actually only created it to avoid uploading excess data from my actual file and extra code. It is therefore not linked to any of my issues. Interestingly, my code somehow worked with the txt as written at first. Anyway, I changed it to csv by now. Additionally, I realized that in Matlab I used slightly different data points. However, all the problems remain. The adjusted csv file is attached, now with the exact same datapoints as were used in Matlab.

The constraints in the fit are actually needed. I have rough estimates where the 3 parameters should be, and anything outside the constraints does not make sense from a physical point of view, especially if the sign is opposite. Also changing the model is not the way to go, as it was derived from the equations corresponding to the physical processes in the experiment.

The standard fit may match very well, however, I need it to be even better. Even with the new csv file, the Matlab fit was far better. Matlab gave: a = 0.0000010000000222047, b = -1.09997779553949E-09, and tau = 657.977318176476. Rsquared was 0.9999992.

Rsquared in Mathematica is 0.9995 when using all the options as I did at first, and a bit lower without them.

I think the fit could reach the same values as Matlab if Mathematica would obey the constraints, which it currently does not for b during the iterations.

Best regards,

Roland

Attachments:
POSTED BY: Roland Käser
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
Attachments
Remove
or Discard

Group Abstract Group Abstract