1
$\begingroup$

I'm trying to use the data:

    In[33]:= data = Transpose[{datat, datax}]

Out[33]= {{1, 83956}, {2, 83959}, {3, 83959}, {4, 83964}, {5, 
  83966}, {6, 83968}, {7, 83970}, {8, 83975}, {9, 83976}, {10, 
  83990}, {11, 84010}, {12, 84011}, {13, 84018}, {14, 84024}, {15, 
  84029}, {16, 84038}, {17, 84044}, {18, 84054}, {19, 84063}, {20, 
  84063}, {21, 84063}, {22, 84063}, {23, 84081}, {24, 84084}, {25, 
  84095}, {26, 84102}, {27, 84103}, {28, 84106}, {29, 84106}, {30, 
  84123}, {31, 84128}, {32, 84146}, {33, 84154}, {34, 84161}, {35, 
  84160}, {36, 84171}, {37, 84177}, {38, 84186}, {39, 84191}, {40, 
  84195}, {41, 84198}, {42, 84209}, {43, 84216}, {44, 84228}, {45, 
  84286}, {46, 84335}, {47, 84378}, {48, 84422}, {49, 84458}, {50, 
  84494}, {51, 84494}, {52, 84553}, {53, 84572}, {54, 84624}, {55, 
  84653}, {56, 84673}, {57, 84701}, {58, 84725}, {59, 84743}, {60, 
  84757}, {61, 84780}, {62, 84785}, {63, 84816}, {64, 84830}, {65, 
  84838}, {66, 84857}, {67, 84871}, {68, 84889}, {69, 84917}, {70, 
  84950}, {71, 84992}, {72, 84992}, {73, 85071}, {74, 85117}, {75, 
  85117}, {76, 85226}, {77, 85246}, {78, 85327}, {79, 85402}, {80, 
  85418}, {81, 85503}, {82, 85622}, {83, 85708}, {84, 85906}, {85, 
  86045}, {86, 86202}, {87, 86381}, {88, 86570}, {89, 86783}, {90, 
  86990}, {91, 87213}, {92, 87489}, {93, 87655}, {94, 87827}, {95, 
  87985}, {96, 88099}, {97, 88206}, {98, 88328}, {99, 88460}, {100, 
  88580}, {101, 88672}, {102, 88793}, {103, 88906}, {104, 
  88958}, {105, 89045}, {106, 89144}, {107, 89214}, {108, 
  89279}, {109, 89375}, {110, 89441}, {111, 89494}, {112, 
  89527}, {113, 89567}, {114, 89616}, {115, 89654}, {116, 
  89695}, {117, 89718}, {118, 89752}, {119, 89784}, {120, 
  89814}, {121, 89836}, {122, 89863}, {123, 89895}, {124, 
  89914}, {125, 89933}, {126, 89953}, {127, 89986}, {128, 
  90008}, {129, 90025}}

To find parameters in the differential equation model that give the best fit:

    pfun = ParametricNDSolveValue[{
   Ii[1.] + S[1.] + R[1.] + F[1.] == 1411780000,
   S'[x] == -R0*F[x]/alpha,
   R'[x] == gamma*Ii[x],
   Ii'[x] == F[x]/alpha + r*F[x] - gamma*Ii[x],
   F'[x] == -F[x]/alpha + R0*F[x]/alpha - r*F[x],
   Ii[1.] == 1900,
   R[1.] == 78504,
   F[1.] == 1900}, R + Ii + F, {x, 0, 150}, {alpha, R0, gamma, r}]

    fit = FindFit[data, 
  pfun[alpha, R0, gamma, r][x], {{alpha, 1.}, R0, gamma, r}, x]

however I got an error code:

 FindFit::nrlnum: The function value {<<1>>} is not a list of real numbers with dimensions {129} at {alpha,R0,gamma,r} = {1.,1.,1.,1.}.

How can I solve this problem? Any help would be greatly appreciated.

$\endgroup$
4
  • 1
    $\begingroup$ One of the most important things when fitting such a complicated function to data is to provide a good initial guess for the parameters. Do you know the approximate values of the parameters (alpha, R0, gamma, r)? At least their order of magnitude? Are they all positive? $\endgroup$
    – Domen
    Jul 28 at 15:30
  • $\begingroup$ Approximately, alpha=14, R0>0, 0<gamma<1, r>0. But I can't find a way to express these using mathematica. $\endgroup$ Jul 28 at 15:38
  • $\begingroup$ I have troubles getting the curve to match your datapoints. Are you sure your ODEs describe the whole range of your data? Can you double check the ODEs and the hard-coded parameters. Take a look at my answer and try it yourself. $\endgroup$
    – Domen
    Jul 28 at 15:59
  • $\begingroup$ It's worth noting that your equations are linear, and so Mathematica can find an exact solution for them without invoking numerical techniques. $\endgroup$ Jul 28 at 16:03
1
$\begingroup$

Apart from @Michael's solution, you can also provide an extra algebraic equation in your system of ODEs:

pfun = ParametricNDSolveValue[{
    Ii[1.] + S[1.] + R[1.] + F[1.] == 1411780000, 
    S'[x] == -R0*F[x]/alpha, R'[x] == gamma*Ii[x], 
    Ii'[x] == F[x]/alpha + r*F[x] - gamma*Ii[x], 
    F'[x] == -F[x]/alpha + R0*F[x]/alpha - r*F[x], Ii[1.] == 1900, 
    R[1.] == 78504, F[1.] == 1900,
    T[x] == R[x] + Ii[x] + F[x]
}, T, {x, 0, 150}, {alpha, R0, gamma, r}];

Now, as I have mentioned in the comments, I have troubles finding the parameters to get a sigmoid-like curve, described by your datapoints. The solution seems to be either exponential-like or hyperbolic-like. Try it yourself:

Manipulate[
 Show[Plot[Evaluate@pfun[alpha, R0, gamma, r]@x, {x, 0, 250},
   PlotRange -> {50000,100000}], ListPlot[data]], {{alpha, 14}, 
  12, 15}, {R0, 0, 2}, {gamma, 0, 1}, {r, 0, 5}]
$\endgroup$
1
$\begingroup$

With pfun from @domen's answer , NonlinearModelFit allows to implement several constraints and evaluates an answer (Method -> "NMinimize" doesn't need starting values)

NonlinearModelFit[data, {pfun[alpha, R0, gamma,r][x], 
alpha > 0, 0 < gamma < 1, R0 > 0, r > 0}, 
{alpha , R0, gamma, r}, x, 
Method -> "NMinimize"]

But the fit is very poor.

$\endgroup$
0
$\begingroup$

The error seems to arise from the fact that you're using ParametricNDSolve to return a result for the anonymous function R + Ii + F, rather than a single anonymous function. This means that when you provide this object with a set of parameters, it returns something of the form

(InterpolatingFunction[data] + InterpolatingFunction[data] +InterpolatingFunction[data])[x]

rather than

InterpolatingFunction[data][x] + InterpolatingFunction[data][x]+InterpolatingFunction[data][x]

as one might expect. There's probably a slick "functional programming" way to fix this using Through or the like, but I'm not enough of a functional programming boffin to know how to do it in this context. (I look forward to seeing if/how other users manage to make this work.)

That said, it's easy enough to recast your ODEs so that they return $f(x) \equiv R(x) + I_i(x) + F(x)$ instead:

eqns = {Ii[1.] + S[1.] + R[1.] + F[1.] == 1411780000, 
  S'[x] == -R0*F[x]/alpha, R'[x] == gamma*Ii[x], 
  Ii'[x] == F[x]/alpha + r*F[x] - gamma*Ii[x], 
  F'[x] == -F[x]/alpha + R0*F[x]/alpha - r*F[x], Ii[1.] == 1900, 
  R[1.] == 78504, F[1.] == 1900};
eqns2 = eqns /. {F[x_] -> f[x] - Ii[x] - R[x], 
  Derivative[n_][F][x_] -> D[f[x] - Ii[x] - R[x], {x, n}]}
pfun2 = ParametricNDSolveValue[eqns2, 
  f, {x, 0, 150}, {alpha, R0, gamma, r}]
fit = FindFit[data, pfun2[alpha, R0, gamma, r][x], {{alpha, 1.}, R0, gamma, r}, x]

This, at least, does not return an error when executed, though it is taking an awfully long time to run on my computer. Providing better initial guesses will no doubt help.

$\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.