Services & ResourcesWolfram ForumsMathGroup Archive

[Date Index] [Thread Index] [Author Index]

Re: SIR epidemic Model

• To: mathgroup at smc.vnet.net
• Subject: [mg36270] Re: [mg36240] SIR epidemic Model
• From: Daniel Lichtblau <danl at wolfram.com>
• Date: Thu, 29 Aug 2002 01:37:54 -0400 (EDT)
• References: <200208280815.EAA18560@smc.vnet.net>
• Sender: owner-wri-mathgroup at wolfram.com

Abdelmajid Khelil wrote:
> 
> Hi everybody,
> 
> I solved the epidemic SIR ODE System
> (<http://library.wolfram.com/webMathematica/MSP/Explore/Biology/Epidemic>
> ) numerically using NDSolve:
> 
> approxsolutions=NDSolve[{
> 
> s´[t]==-a s[t] i[t],
> i´[t]==a s[t] i[t] - b i[t],
> r´[t]==b i[t],
> 
> s[0]==700, i[0]==1, r[0]==0}, {i[t], s[t], r[t]}, {t,0,20}];
> 
> and this for fixed values of a and b.
> 
> Now I want to fit the solution for variable a and b to a list of given
> points (fitting using the least squares method and the mathematica fonction
> findminimum). My questions are
> 1- Is it possible to solve the ODE (with mathematica) for not fixed values
> of a and b (so to get a parametric solution that depends on a and b)?
> 2- Is it possible in mathematica to use the fonction findminimum with the
> (not explicit) solutions of the ODE, if not are there some other
> alternatives?
> 
> In advance, thank you very much
> 
> Majid


There was a discussion on this topic three years ago. The URL below will
take you to a reply in that thread that includes remarks and code from
Carl Woll and myself. I'll also show some code for your specific
example, in order to indicate potential problems.

http://library.wolfram.com/mathgroup/archive/1999/Oct/msg00104.html

Below is some mildly clumsy code tailored to your example. It will do
the DE solving given parameter values, form an objective function
(sum-of-squares type) given actual data and solutions for specific
parameters, and to do the minimization. I did everything at machine
precision but this is not necessarily a good idea. Also I do no error
checking e.g. for cases where the numeric DE solver fails to solve over
the full range. A related issue is that, depending on roughly where the
correct parameter values live, and how close you are in your initial
guess, this might not be at all a well-behaved problem.

machsoln[a_?NumericQ,b_?NumericQ] := First[NDSolve[{
	D[ss[t],t]==-a*ss[t]*ii[t], D[ii[t],t]==a*ss[t]*ii[t]-b*ii[t],
	  D[rr[t],t]==b*ii[t], ss[0]==700, ii[0]==1, rr[0]==0},
	  {ii[t], ss[t], rr[t]}, {t,0,20}]];

machsol = machsoln[1/100,1];

We'll generate some data for parameter values a=1/100, b=1.

machdata = Table[{t,ii[t],ss[t],rr[t]} /. machsol, {t,0,20,5}]

For the objective function, given numeric values for {a,b}, we'll solve
the system, evaluate at the same 't' values as are used for our data,
subtract resulting vectors from corresponding data vectors, take norm
squares, and sum.

objfunc[data_,{a_?NumericQ,b_?NumericQ}] := Module[
	{iii, sss, rrr, newdata, vals},
	{iii,sss,rrr} = {ii[t],ss[t],rr[t]} /. machsoln[a,b];
	newdata = Map[{iii,sss,rrr} /. t-># &, Map[First,data]];
	vals = Map[#.#&, Map[Rest,data]-newdata];
	Apply[Plus, vals]
	]

For example, we will get 0 if we use {a,b} = {.01,1}, and something
nonzero if we increase both modestly, say by 10%.

In[85]:= objfunc[machdata, {.01,1.}]
Out[85]= 0.

In[86]:= objfunc[machdata, {.011,1.1}]
Out[86]= 72.5289

With the above objective function we can call FindMinimum. For this we
will use two starting values for each parameter so that FIndMinimum can
use a secant method. An alternative to this, wherein one passes along a
gradient "black box" evaluator, is presented by Carl Woll in his note at
the URL above.

findParams[data_, inits:{a1_,b1_}] := Module[
	{vals,a,b},
	FindMinimum[objfunc[data,{a,b}], {a,a1,a1+.1*a1}, {b,b1,b1+.1*b1}]
	]

For example, if I give {a,b} = {.011,1.1} as initial values, I get a
result reasonably fast and quite accurate, albeit with many warnings
indicated. They basically indicate 

In[87]:= findParams[data, {.011,1.1}]

InterpolatingFunction::dmval: 
   Input value {20} lies outside the range of data in the interpolating
     function. Extrapolation will be used.

InterpolatingFunction::dmval: 
   Input value {20} lies outside the range of data in the interpolating
     function. Extrapolation will be used.

InterpolatingFunction::dmval: 
   Input value {20} lies outside the range of data in the interpolating
     function. Extrapolation will be used.

General::stop: Further output of InterpolatingFunction::dmval
     will be suppressed during this calculation.

                    -11
Out[87]= {8.46495 10   , {a$4965 -> 0.01, b$4965 -> 1.}}

It still does reasonably well when I perturb by 50% from the "correct"
parameter values. This time the warning message indicates that we might
be closer to trouble.

In[89]:= findParams[data, {.015,1.5}]

FindMinimum::fmcv: 
   Failed to converge to the requested accuracy or precision within 30
     iterations.

                    -10
Out[89]= {6.48159 10   , {a$8747 -> 0.01, b$8747 -> 1.}}

You might want to play with various options in FindMinimum such as
MaxIterations. You might also need to use NDSolve options e.g.
WorkingPrecision in order to improve the quality of the result.


Daniel Lichtblau
Wolfram Research

• References:
  • SIR epidemic Model
    • From: "Abdelmajid Khelil" <khelil@informatik.uni-stuttgart.de>