# Problem with NonlinearModelFit in the case of ODE

I am asking for your help for fitting a nonlinear ODE to data which are unknown up to a multiplicative constant k. So the constant k appears both in the data, and in the model, where it intervenes as initial condition. There is no problem for having k appear only in the initial condition -- see the attached documentation example, where we added k

v1=Table[i,{i,0,32}];
v2={107, 135, 612, 195, 626, 619, 491, 1164, 1137, 511, 1036, 1144, 2650, 3162, 6074, 6693,
8253, 6639, 6148, 4345, 3141, 1958, 1130, 484, 356, 296, 195, 121, 208, 101,
67, 128, 2};(*initial values of infections*)

data=Inner[List,v1,k v2, List]

model[\[Gamma]_?NumberQ,\[Beta]_?NumberQ, s0_?NumberQ, k_?NumberQ]:=model[\[Gamma],\[Beta],s0,k]=NDSolveValue[
{s'[t]==-\[Beta] i[t] *s[t],i'[t]==-(\[Gamma]) i[t]+\[Beta] i[t] *s[t],s[0]==s0,i[0]==  data[[1,2]]},{s,i},{t,0,33}][[1]]

nlm=NonlinearModelFit[ data,model[\[Gamma],\[Beta],s0,k][t],{\[Gamma], \[Beta],s0,k},t]

par=nlm["BestFitParameters"]

{tmin,tmax}=MinMax[Inner[List,v1, v2, List][[All,1]]];
Show[ListPlot[Map[{#[[1]], #[[2]]}&/. par,Inner[List,v1, v2, List]]],
Plot[nlm[t],{t,tmin,tmax},PlotStyle->Orange],Frame->True,PlotRange->All]

But, when I add k also in the data ... the syntax is accepted, nlm still works, but the fit is very bad. Any help is gratefully received. Thanks.

• What are the t-values in data? Jan 24, 2023 at 14:32
• Dear @UlrichNeumann , you will find above the x-values corresponding to the data within the data coordinates. Could you please, clarify why you would require these values, and what is their relevance to this problem? Jan 24, 2023 at 17:27
• Usually NonlinearModelFitexpects an array. If data is onedimensional the indices of data are assumed to be the x-values! Jan 25, 2023 at 10:04
• In your latest edit of the question, k has disappeared in the definition of the model. Also, it appears that when the response variable has the form {107 k, 135 k, 612 k, 195 k,...} that k is treated as a separate entity from the k in the parameter list. I'm not aware that the first argument in NonlinearModelFit corresponding to the data can be anything but numbers and no symbols. Do you have a reference for why that's not true?
– JimB
Jan 27, 2023 at 17:29

• Since we want to fit the second function i[t], so we need to use [[2]] instead of [[1]].

• Method -> "LevenbergMarquardt" seems to work for the original data( without k though many warning message)

Clear["Global`*"];
data = {{0, 107}, {1, 135}, {2, 612}, {3, 195}, {4,
626}, {5, 619}, {6, 491}, {7, 1164}, {8, 1137}, {9, 511}, {10,
1036}, {11, 1144}, {12, 2650}, {13, 3162}, {14, 6074}, {15,
6693}, {16, 8253}, {17, 6639}, {18, 6148}, {19, 4345}, {20,
3141}, {21, 1958}, {22, 1130}, {23, 484}, {24, 356}, {25,
296}, {26, 195}, {27, 121}, {28, 208}, {29, 101}, {30, 67}, {31,
128}, {32, 20}};(*initial values of infections*)
model[γ_?NumericQ, β_?NumericQ, s0_?NumericQ,
k_?NumericQ] :=
NDSolveValue[{s'[t] == -β i[t]*s[t],
i'[t] == -(γ) i[t] + β i[t]*s[t], s[0] == s0,
i[0] == k*data[[1, 2]]}, {s, i}, {t, 0, 33}][[2]];
nlm = NonlinearModelFit[data,
model[γ, β, s0, k][t], {γ, β, s0, k}, t,
Method -> "LevenbergMarquardt"]
par = nlm["BestFitParameters"]
{xmin, xmax} = MinMax[data[[All, 1]]];
Show[ListPlot[data],
Plot[nlm[t], {t, xmin, xmax}, PlotStyle -> Orange], Frame -> True,
PlotRange -> All]

• Thanks for your help. Jan 25, 2023 at 19:34

Your model which I'm led to believe is

p * model[γ, β, s0, k][t]

has numerical instability issues. (It would not be too extreme to say that your model with the inclusion of the proportionality parameter (which I've included as p) is close to being overparameterized. (Also, you mention k as the proportionality parameter but you've already got a k elsewhere in the model. Are you wanting those two values to be identical?)

Part of the instability is caused by the large orders of magnitude differences in the values of the parameters. Functions that use an iterative process to find a minimum or maximum don't always work so well when the parameters differ by large orders of magnitude. Further your code uses the starting values of 1 for all parameters. Using the estimates from @cvgmt 's answer as the starting values gets a slightly better fit in that the estimate of the standard error of estimate is reduced slightly from 451.276 to 448.935:

nlm["EstimatedVariance"]^0.5
(* 451.276 *)

nlm1 = NonlinearModelFit[data,
model[γ, β, s0, k][
t], {{γ, 815/1000}, {β, 266/10000000}, {s0,
579561/10}, {k, 366/100000}}, t, Method -> "LevenbergMarquardt",
WorkingPrecision -> 30] // Quiet;
nlm1["EstimatedVariance"]^0.5
(* 427.609 *)

Now adding in a proportionality parameter actually increases the resulting standard error of estimate:

nlm2 = NonlinearModelFit[data,
p model[γ, β, s0, k][t], {{γ,
815/1000}, {β, 266/10000000}, {s0, 579561/10}, {k,
366/100000}, {p, 1}}, t, Method -> "LevenbergMarquardt",
WorkingPrecision -> 20] // Quiet
nlm2["EstimatedVariance"]^0.5
(* 432.65 *)

This increase with more parameters is another indication of numerical instability.

Further, the more complex model increases the $$AIC_c$$ value:

#["AICc"]^0.5 & /@ {nlm, nlm1, nlm2}
(* {22.4782, 22.399, 22.461} *)

I recommend consulting a statistician as thinking that the regression is a multiple of a solution to a differential equation seems a bit arbitrary. And if being arbitrary is OK, then a nonparametric regression would likely provide a better fit.

One indicator of an overparameterized model is that the correlations of the parameter estimators are near +1 or -1. That is the case here. Consider the correlation matrix of @cvgmt 's model:

nlm["CorrelationMatrix"] // MatrixForm

(And this is before including a proportionality parameter which will make the results more unstable.) This doesn't mean that the model is wrong but given the available data, the parameter estimators are highly correlated and have large individual standard errors. This means that very different sets of parameters result in nearly the same predictions.

Even putting in the results of that model as starting values results in wildly different parameter values but has nearly the same predictions.

This is not an issue or error with NonlinearModelFit or NDSolveValue. And using a higher value for WorkingPrecision can't fix it.

In short, the predictions are fine but one can't pin down the individual parameters very well.

Sometimes (many times?) an iterative algorithm stops too soon. (Again, my experience is that this happens more often when parameters to be estimated differ by many orders of magnitude.) One approach is to restart the process again with the estimates of the previous result and repeat that until there is stabilization of the estimates: lather, rinse, repeat. Here I start with @cvgmt 's code and repeat a few times.

nlm = NonlinearModelFit[data, model[γ, β, s0, k][t], {γ, β, s0, k}, t,
Method -> "LevenbergMarquardt"] // Quiet;

{γ2, β2, s02, k2} = {γ, β, s0, k} /. nlm["BestFitParameters"];
nlm2 = NonlinearModelFit[data, model[γ, β, s0, k][t],
{{γ, γ2}, {β, β2}, {s0, s02}, {k, k2}}, t,
Method -> "LevenbergMarquardt"] // Quiet;

{γ3, β3, s03, k3} = {γ, β, s0, k} /. nlm2["BestFitParameters"];
nlm3 = NonlinearModelFit[data, model[γ, β, s0, k][t],
{{γ, γ3}, {β, β3}, {s0, s03}, {k, k3}}, t,
Method -> "LevenbergMarquardt"] // Quiet;

{γ4, β4, s04, k4} = {γ, β, s0, k} /. nlm2["BestFitParameters"];
nlm4 = NonlinearModelFit[data, model[γ, β, s0, k][t],
{{γ, γ4}, {β, β4}, {s0, s04}, {k, k4}}, t,
Method -> "LevenbergMarquardt"] // Quiet;

The parameter estimates vary considerably although the overall fit looks pretty much the same:

#["ParameterTable"] & /@ {nlm, nlm2, nlm3, nlm4} // TableForm

{xmin, xmax} = MinMax[data[[All, 1]]];
Show[ListPlot[data], Plot[{nlm[t], nlm2[t], nlm3[t]}, {t, xmin, xmax},
PlotStyle -> {{Thickness[0.02], Orange}, {Thickness[0.01], Blue},       Yellow},
PlotLegends -> "Expressions"], Frame -> True, PlotRange -> All]

(Note that the fit for nlm4 is not included as it is the same as the fit for nlm3.)

The estimated standard errors of estimate are smallest with nlm3 and nlm4:

#["EstimatedVariance"]^0.5 & /@ {nlm, nlm2, nlm3, nlm4} // TableForm

And the parameter correlation matrices look less problematic with the final results:

(# // MatrixForm &) /@ (#["CorrelationMatrix"] & /@ {nlm, nlm2, nlm3, nlm4})

• Thanks, @cvgmt and JimB. This was of great help. It seems the weak point in our attempt is that 33 weekly data are not enough to provide reliable ODE modelling for an epidemic. Maybe it will make more sense with daily data Jan 30, 2023 at 11:15
• @florin: Please note that my comments are about the statistical aspects as you and Rim way outclass me in the mathematics. Determining the adequacy of the amount and type of data and model can be best done with simulations on known values of the parameters. And because you have "counts", the error is probably best modeled as a Poisson or negative binomial (which NonlinearModelFit does not allow). Also, given the fit from cvgmt, I think there is a lack of fit especially at the beginning of the sequence. Here, too, simulations will help you determine that.
– JimB
Jan 30, 2023 at 15:09
• Dear @JimB, thanks for your edited answer, which helped display the best-fitted parameters much better. It seems more evident with this method that the estimated parameters start looking similar once we reduce the data capacity (for instance, to 25, in this case, nlm2, nlm3 and nlm4 look similar). Jan 31, 2023 at 11:53
• Dear @JimB, we are indeed mathematicians (now in epidemiology), and this was our first attempt to confront data. Your help was invaluable. We would like to work further in this direction under your guidance, if you might have time. Can we contact you by email? Feb 2, 2023 at 14:20
• @florin Certainly. If the Pau address good for you?
– JimB
Feb 2, 2023 at 15:51

data = {{6.47, 3.65}, {7.43,
3.45}, {3.9, -2.94}, {4.8, -1.29}, {2.48, -0.35}, {6.32,
3.16}, {2.59, -1.19}, {9.13, -2.}, {3.81, -3.04}, {3.33, -2.68}};

model[a_?NumberQ, b_?NumberQ, c_?NumberQ, k_?NumberQ] :=
Module[{y, x},First[y /. NDSolve[{y''[x] + a y[x] == 0, y[0] ==
b*k*data[[1,1]],y'[0] == c}, y, {x, 0, 10}]]]

nlm = NonlinearModelFit[data, model[a, b, c, k][x], {a, b, c, k}, x,Method -> "Gradient"]