# Parameter Estimation on ODE Using NMinimize (Least Squares Method)

I'm trying to estimate the parameters in this equation using this code:

data1 = {0.006,0.0182,0.0191,0.0154,0.0323,0.0199,0.0339,0.016,0.0309,0.0291,0.0279,0.0308,0.0353,0.0235,0.0223,0.0201,0.0311,
0.0405,0.0284,0.0252,0.0186,0.0218,0.0176,0.033,0.0467,0.0535,0.0461,0.0451,0.0281,0.0309,0.0399,0.0462,0.0439,0.0413,0.049,
0.0321,0.0439,0.0529,0.0584,0.0461,0.0491,0.0398,0.0341,0.0382,0.0509,0.052,0.0476,0.0411,0.0337,0.0261,0.0335,0.0391,0.0359,
0.0334,0.0253,0.026,0.0112,0.0126,0.024,0.0232,0.0207,0.0179,0.0142,0.0044,0.0153,0.0124,0.0171,0.0252,0.0183,0.0142,0.0122,0.0165,
0.0101,0.0195,0.011,0.0092,0.009,0.0075,0.0097,0.0068,0.0077,0.0088,0.0107,0.004,0.0046,0.0051,0.0044,0.0063,0.0074,0.0052,0.0074,
0.0019,0.0075,0.0047,0.0051,0.005,0.0016,0.0041,0.0028,0.0075,0.0041,0.0051,0.006,0.0095,0.003,0.0024,0.006,0.0031,0.0044,0.0039,
0.005,0.0033,0.0025,0.0024,0.0038,0.0028,0.0041,0.0034,0.0027,0.0032,0.0035,0.005};
T = Length[data1] - 1;

(*Parameters*)
Lambda=0.5;
Alpha=0.5;
Kappa=0.00398;
d=0.0047876;
r=0.09871;
tau=7;

SEIRVB[Lambda_, Alpha_, Kappa_, d_, r_,tau_] := NDSolve[
{S'[t]==Lambda - Beta1 S[t] E1[t] - (Beta3 + Alpha) S[t],
E1'[t]==Beta4 V[t] E1[t] + Beta6 B[t] E1[t] + Beta1 S[t] E1[t] - Beta2 E1[t] I1[t] - (Kappa + Alpha) E1[t],
I1'[t]==Beta5 V[t-tau] I1[t] + Beta7 B[t-tau] I1[t] + Beta2 E1[t] I1[t] - (d + Alpha + r) I1[t],
R'[t]==Beta9 V[t] + Beta8 B[t] + r I1[t] + Kappa E1[t] - Alpha R[t],
V'[t]==Beta3 S[t] - (Beta9 + Alpha + Beta10) V[t] - Beta4 V[t] E1[t] - Beta5 V[t] I1[t - tau],
B'[t]==Beta10 V[t] - Beta6 B[t] E1[t] - Beta7 B[t-tau] I1[t] - (Beta8 + Alpha) B[t],
S[0]==0.3,
E1[0]==0.5,
I1[0]==0.06,
R[0]==0,
V[0]==0,
B[0]==0},
{S, E1, I1, R, V, B},
{t,0,T},
Method -> "ImplicitRungeKutta"];

leastSquaresObjective[Lambda_, Alpha_, Kappa_, d_, r_, tau_,infectedData_] :=
Module[{sol, t, infectedEstimate, infectedResiduals},
sol = SEIRVB[Lambda, Alpha, Kappa, d, r, tau];
infectedEstimate = I[t] /. sol;
infectedResiduals = infectedData - infectedEstimate;
Total[infectedResiduals^2]
]

result = Minimize[
{leastSquaresObjective[Lambda, Alpha, Kappa, d, r, tau, data1],
Beta1>0,Beta2>0,Beta3>0,Beta4>0,Beta5>0,Beta6>0,Beta7>0,Beta8>0,Beta9>0,Beta10>0},
{Beta1, Beta2, Beta3, Beta4, Beta5, Beta6, Beta7, Beta8, Beta9, Beta10}
]


However I always got an error with this:

General:: ReplaceAll::reps: {NDSolve[{(S^[Prime])[t]==0.5 -(0.5 +[Beta]3) S[t]-[Beta]1 E1[t] S[t],(E1^[Prime])[t]==-0.50398 E1[t]+[Beta]6 B[t] E1[t]-[Beta]2 E1[t] I1[t]+[Beta]1 E1[t] S[t]+[Beta]4 E1[t] V[t],(I1^[Prime])[t]==-0.603498 I1[t]+[Beta]7 B[Plus[<<2>>]] I1[t]+[Beta]2 E1[t] I1[t]+[Beta]5 I1[t] V[Plus[<<2>>]],(R^[Prime])[t]==[Beta]8 B[t]+0.00398 E1[t]+0.09871 I1[t]-0.5 R[t]+[Beta]9 V[t],(V^[Prime])[t]==[Beta]3 S[t]-(0.5 +[Beta]10+[Beta]9) V[t]-[Beta]4 E1[t] V[t]-[Beta]5 I1[Plus[<<2>>]] V[t],(B^[Prime])[t]==-((0.5 +[Beta]8) B[t])-[Beta]6 B[t] E1[t]-[Beta]7 B[Plus[<<2>>]] I1[t]+[Beta]10 V[t],S[0]==0.3,E1[0]==0.5,I1[0]==0.06,R[0]==0,<<2>>},{S,E1,I1,R,V,B},{t,0,121},Method->ImplicitRungeKutta]} is neither a list of replacement rules nor a valid dispatch table, and so cannot be used for replacing.

Along with this error in NMinimize:

NMinimize::nnum: The function value {(0.0016 -I[t$$117055])^2+(0.0019 -I[t$$117055])^2+2 (0.0024 -I[t$$117055])^2+(0.0025 -I[t$$117055])^2+(0.0027 -I[t$$117055])^2+2 (0.0028 -I[t$$117055])^2+(0.003 -I[t$$117055])^2+(0.0031 -I[t$$117055])^2+(0.0032 -I[t$$117055])^2+(0.0033 -I[t$$117055])^2+<<92>>} is not a number at {[Beta]1,[Beta]10,[Beta]2,[Beta]3,[Beta]4,[Beta]5,[Beta]6,[Beta]7,[Beta]8,[Beta]9} = {1.23676,1.7426,1.81554,1.02898,0.705467,1.51602,1.23152,1.60226,1.40791,0.936109}.

And after that, even though I replaced the values of kappa, r and the data itself, the end result is always the NMinimize error one.

Does anyone have an ide on how to fix this? Thanks in advance!

• The code for NDSolve has syntax errors, please fix them. Commented Jun 18, 2023 at 10:31
• Fixed it already, apologies Commented Jun 18, 2023 at 14:38
• The endpoint of the integration T is not defined. Please test your own code. Commented Jun 18, 2023 at 16:01
• @Malanie Could you give some explanation to your model? Commented Jun 18, 2023 at 16:13
• I've given a provisional answer below. I would like to point out that parametric fitting of discrete data with ODE's is the power zone of Pontryagin's Minimum Principle, where the performance index is sum of squares of observed-expected. Commented Jun 23, 2023 at 23:28

1. You will not be able to use I[t] as capital I is reserved for the imaginary number. This may be I1 instead.
2. You will have to pass Beta's to the SEIRVB and leastSquaresObjective, likely as Beta1_?NumericQ, Beta2_?NumericQ, ... to prevent symbolic Beta1, Beta2's, etc. NDSolve needs to evaluate to numerical values everywhere.
3. It is necessary to evaluate the interpolated function at each time in {t,0,T}.
4. Please note that t - tau is a negative number for t < tau. This may not be an issue, but Mathematica will issue warnings.
5. I have utilized the Automatic method rather than "ImplicitRungeKutta" as the latter seems to take a while.

Here is my attempt to edit your code:

data1 = {0.006, 0.0182, 0.0191, 0.0154, 0.0323, 0.0199, 0.0339,
0.016, 0.0309, 0.0291, 0.0279, 0.0308, 0.0353, 0.0235,
0.0223, 0.0201, 0.0311, 0.0405, 0.0284, 0.0252, 0.0186,
0.0218, 0.0176, 0.033, 0.0467, 0.0535, 0.0461, 0.0451,
0.0281, 0.0309, 0.0399, 0.0462, 0.0439, 0.0413, 0.049,
0.0321, 0.0439, 0.0529, 0.0584, 0.0461, 0.0491, 0.0398,
0.0341, 0.0382, 0.0509, 0.052, 0.0476, 0.0411, 0.0337,
0.0261, 0.0335, 0.0391, 0.0359, 0.0334, 0.0253, 0.026,
0.0112, 0.0126, 0.024, 0.0232, 0.0207, 0.0179, 0.0142,
0.0044, 0.0153, 0.0124, 0.0171, 0.0252, 0.0183, 0.0142,
0.0122, 0.0165, 0.0101, 0.0195, 0.011, 0.0092, 0.009,
0.0075, 0.0097, 0.0068, 0.0077, 0.0088, 0.0107, 0.004,
0.0046, 0.0051, 0.0044, 0.0063, 0.0074, 0.0052, 0.0074,
0.0019, 0.0075, 0.0047, 0.0051, 0.005, 0.0016, 0.0041,
0.0028, 0.0075, 0.0041, 0.0051, 0.006, 0.0095, 0.003,
0.0024, 0.006, 0.0031, 0.0044, 0.0039, 0.005, 0.0033,
0.0025, 0.0024, 0.0038, 0.0028, 0.0041, 0.0034, 0.0027,
0.0032, 0.0035, 0.005};
T = Length[data1] - 1;

(*Parameters*)
Lambda = 0.5;
Alpha = 0.5;
Kappa = 0.00398;
d = 0.0047876;
r = 0.09871;
tau = 7;

SEIRVB[Lambda_, Alpha_, Kappa_, d_, r_, tau_, Beta1_?NumericQ,
Beta2_?NumericQ, Beta3_?NumericQ, Beta4_?NumericQ, Beta5_?NumericQ,
Beta6_?NumericQ, Beta7_?NumericQ, Beta8_?NumericQ,
Beta9_?NumericQ, Beta10_?NumericQ] :=
NDSolve[{S'[t] == Lambda - Beta1 S[t] E1[t] - (Beta3 + Alpha) S[t],
E1'[t] ==
Beta4 V[t] E1[t] + Beta6 B[t] E1[t] + Beta1 S[t] E1[t] -
Beta2 E1[t] I1[t] - (Kappa + Alpha) E1[t],
I1'[t] ==
Beta5 V[t - tau] I1[t] + Beta7 B[t - tau] I1[t] +
Beta2 E1[t] I1[t] - (d + Alpha + r) I1[t],
R'[t] ==
Beta9 V[t] + Beta8 B[t] + r I1[t] + Kappa E1[t] - Alpha R[t],
V'[t] ==
Beta3 S[t] - (Beta9 + Alpha + Beta10) V[t] - Beta4 V[t] E1[t] -
Beta5 V[t] I1[t - tau],
B'[t] ==
Beta10 V[t] - Beta6 B[t] E1[t] -
Beta7 B[t - tau] I1[t] - (Beta8 + Alpha) B[t], S[0] == 0.3,
E1[0] == 0.5, I1[0] == 0.06, R[0] == 0, V[0] == 0, B[0] == 0}, {S,
E1, I1, R, V, B}, {t, 0, T}, Method -> Automatic (*"ImplicitRungeKutta"*)];

leastSquaresObjective[Lambda_, Alpha_, Kappa_, d_, r_, tau_,
infectedData_, Beta1_?NumericQ, Beta2_?NumericQ, Beta3_?NumericQ,
Beta4_?NumericQ, Beta5_?NumericQ, Beta6_?NumericQ, Beta7_?NumericQ,
Beta8_?NumericQ, Beta9_?NumericQ, Beta10_?NumericQ] :=
Module[{sol, t, infectedEstimate, infectedResiduals},
sol = SEIRVB[Lambda, Alpha, Kappa, d, r, tau, Beta1, Beta2, Beta3,
Beta4, Beta5, Beta6, Beta7, Beta8, Beta9, Beta10];
infectedEstimate = First[I1 /. sol];
infectedResiduals =
infectedData - Table[infectedEstimate[t], {t, 0, T}];
Total[infectedResiduals^2]];

result =
NMinimize[{leastSquaresObjective[Lambda, Alpha, Kappa, d, r, tau,
data1, Beta1, Beta2, Beta3, Beta4, Beta5, Beta6, Beta7, Beta8,
Beta9, Beta10], Beta1 > 0, Beta2 > 0, Beta3 > 0, Beta4 > 0,
Beta5 > 0, Beta6 > 0, Beta7 > 0, Beta8 > 0, Beta9 > 0,
Beta10 > 0}, {Beta1, Beta2, Beta3, Beta4, Beta5, Beta6, Beta7,
Beta8, Beta9, Beta10}]


The result:

{0.0261644, {Beta1 -> 3.51855, Beta2 -> 0.721769, Beta3 -> 0.,
Beta4 -> 2.2859, Beta5 -> 2.30187, Beta6 -> 1.07471,
Beta7 -> 0.00360494, Beta8 -> 0.125174, Beta9 -> 2.95638,
Beta10 -> 1.93584}}


It is also use to make plots:

Print[ListPlot[{infectedData,
Table[infectedEstimate[t], {t, 0, T}]}]];


of the experimental data vs the simulated data inside the leastSquaresObjective function, and so at the fitted parameters:

I doubt this is a reasonable solution, and you might wish to try other starting values to check for a better solution.

Scanning through all the plots for the NMinimize steps, I did see:

which has perhaps closer features to the experimental data, but still not quite right.

It is also advisable to double-check the formulae for the differential equation and all of the other steps to make sure there are no errors. This is merely a starting point that will not error out right away.

You need to use ParametricNDSolve with Beta1, Beta2,...,Beta10 as parameters. NDSolve won't work with non-numerical parameters:

SEIRVB[Lambda_, Alpha_, Kappa_, d_, r_, tau_] :=
ParametricNDSolve[{S'[t] ==
Lambda - Beta1 S[t] E1[t] - (Beta3 + Alpha) S[t],
E1'[t] ==
Beta4 V[t] E1[t] + Beta6 B[t] E1[t] + Beta1 S[t] E1[t] -
Beta2 E1[t] I1[t] - (Kappa + Alpha) E1[t],
I1'[t] ==
Beta5 V[t - tau] I1[t] + Beta7 B[t - tau] I1[t] +
Beta2 E1[t] I1[t] - (d + Alpha + r) I1[t],
R'[t] ==
Beta9 V[t] + Beta8 B[t] + r I1[t] + Kappa E1[t] - Alpha R[t],
V'[t] ==
Beta3 S[t] - (Beta9 + Alpha + Beta10) V[t] - Beta4 V[t] E1[t] -
Beta5 V[t] I1[t - tau],
B'[t] ==
Beta10 V[t] - Beta6 B[t] E1[t] -
Beta7 B[t - tau] I1[t] - (Beta8 + Alpha) B[t], S[0] == 0.3,
E1[0] == 0.5, I1[0] == 0.06, R[0] == 0, V[0] == 0, B[0] == 0}, {S,
E1, I1, R, V, B}, {t, 0, T}, {Beta1, Beta2, Beta3, Beta4, Beta5,
Beta6, Beta7, Beta8, Beta9, Beta10},
Method -> "ImplicitRungeKutta"]


Then we can call SEIRVB and collect the parametric solution functions (I call the solution for I1 infectedEstimate because I believe that's what you meant in the part of the question infectedEstimate = I[t] /. sol but correct me if I misinterpreted this):

{Ssoln, E1soln, infectedEstimate, Rsoln, Vsoln,
Bsoln} = {S, E1, I1, R, V, B} /.
SEIRVB[Lambda, Alpha, Kappa, d, r, tau];



Now I'm assuming the time coordinates for the data in data1 is {0,1,...,T} since you defined T=Length[data1] -1. So we can define our objective function as:

squareDiff[parList_] :=
Module[{f = infectedEstimate[parList /. List -> Sequence]},

fTab = f[Range[0, T]];

Norm[data1 - fTab]
]


This is the part I cannot get to work though, maybe someone else will have an idea why

varList = Array[beta, 10];
NMinimize[{squareDiff[varList], varList \[VectorGreater] 0}, varList]



This gives the error message (very similar to your original NMinimize error):

NMinimize::nnum: The function value {0.47362,0.591819,0.594828,0.506271,0.381692,0.267308,0.193566,0.173022,0.182947,0.199315,<<112>>} is not a number at {beta[1],beta[2],beta[3],beta[4],beta[5],beta[6],beta[7],beta[8],beta[9],beta[10]} = {1.23676,1.7426,1.81554,1.02898,0.705467,1.51602,1.23152,1.60226,1.40791,0.936109}.

However, plugging {1.23676,1.7426,1.81554,1.02898,0.705467,1.51602,1.23152,1.60226,1.40791,0.936109}

into squareDiff does indeed give a number:

squareDiff[{1.23676, 1.7426, 1.81554, 1.02898, 0.705467, 1.51602,
1.23152, 1.60226, 1.40791, 0.936109}]

(*0.289701*)


Maybe someone else know's what is causing this error though. I suspect I defined squareDiff in a non-optimal way that could be causing NMinimize to think it is not numerically valued at certain input values.