Message Boards

GROUPS:

Robert Rimmer

[NB] Logistic Model for Quarantine Controlled Epidemics

Robert Rimmer

Posted 2 years ago

10396 Views

6 Replies

8 Total Likes

Follow this post

MODERATOR NOTE: coronavirus resources & updates: https://wolfr.am/coronavirus

A logistic growth model tracks the Coronavirus COVID-19 epidemic, using the data provided by Johns Hopkins CSSE and uploaded to GitHub. The data are also available from the Wolfram Data Repository. The model is very simple and follows the function proposed by Pierre François Verhulst 180 years ago. The original concept behind the function was to model population growth, which should be limited by lack of resources and space. Epidemic dynamics, expressed as a cumulative number of cases or deaths, can use the same model, when the primary method of control is quarantine as in the case of a novel viral infection.

There are only three parameters to the function. The parameter, k, is the growth rate factor; this parameter is determines the rate of infection spread and how long the epidemic will last. The parameter, L, is the limit of the number of infected people or the subset of people who die. The t0 parameter is the time of the peak of the epidemic on an arbitrary time scale.

MoreinformationandhistoryonthelogisticfunctionfromWikipedia.

Readers, who do not wish to get into the mathematical weeds immediately, can look at the structure of the function and use the Manipulate cells at the end of the next section, then read “Why the Model Works” and jump to the last section which analyzes some actual data.

This notebook does not link to the databases, but contains all the Mathematica functions used to analyze the data. The goal is to explore the logistic function in detail, demonstrate that it is a good model for quarantine limited epidemics, and show some examples of how the model can be used to track an epidemic.

The Derivation of the Logistic Function

Logistic growth is defined by a differential equation. Let f(t) be the cumulative count of infected cases at time, t.

f(t)

t

kf(t)1-

f(t)

where k is a continuous growth rate and L is the limit to growth such that f(t) cannot exceed L. When f(t) is small relative to L, the term,

f(t)

is near 1. Then the initial growth is nearly proportional to f(t) with exponential growth rate, k. When f(t) reaches L, the growth rate is zero and growth stops. If L and k are positive, then the growth rate is always positive.

The DSolve function solves this simple differential equation.

DSolve[{f'[t]kf[t](1-f[t]/L),f[t0]L/2},f[t],t]//FullSimplify

Out[]=

f[t]



kt0





This can be simplified:

Exp[-kt]NumeratorDenominator



kt0



//Expand#1/#2&@@%//Simplify

FullSimplify



kt0





k(-t+t0)





This proves the solution using the differential equation:

D

k(-t+t0)



,tk

k(-t+t0)



k(-t+t0)



L//FullSimplify

The boundary condition, f(t0)  L/2, sets t0 to be the inflection point where growth rate is maximum, i.e. the time of the peak of the epidemic.

SolveD

k(-t+t0)



,{t,2}0,t,Reals

The manipulate cell shows how the parameters affect the logistic function. L simply changes the scale on the y-axis, t0 changes the location of the center of the function on the x-axis. The logistic function is symmetric horizontally and vertically about the point, {t0, L/2}, which is set as the axis origin on the logistic function graph. k is the rate parameter which determines the duration of the epidemic.

ManipulateGraphicsColumnPlot

-k(t-t0)



,{t,0,60},AxesOrigin{t0,L/2},AxesLabel{"Time","Cumulative Cases"},PlotRangeAll,PlotLabel"Logistic Function",Plot

k(t+t0)



(



kt0



)

,{t,0,60},PlotStyleLighter@Orange,PlotRangeAll,AxesLabel{"Time","Cases per Unit Time"},PlotLabel"First Derivative Logistic Function",{{L,1000},500,100000},{{k,0.3},0.1,1},{{t0,30},20,40},ContinuousActionTrue,SaveDefinitionsTrue

Out[]=

The denominator of the logistic function, in orange, approaches the limit, 1, as time increases. As the denominator approaches 1, the function approaches the limit, L. The growth rate parameter, k, controls how quickly the limit is reached. For the epidemiologic data the time scale will be in days.

Manipulate[Plot[Evaluate[

-k(t-t0)



,1+

-k(t-t0)



/.{t00,k0.5}],{t,-1,50},PlotLegends"Expressions",GridLines{None,{{1,{Black,Thick}}}},PlotLabel"Logistic Denominator",AxesLabel{"Time",None}],{{k,0.5},0.01,1},SaveDefinitionsTrue,ContinuousActionTrue]

Out[]=

Assuming the epidemic begins when the logistic function equals one, the duration of the epidemic can be calculated.

sol=RefineSolve

-k(t-t0)



1,t,Reals,L>1//PowerExpand//Simplify

Out[]=

tt0-

Log[-1+L]



By symmetry, where t0 is the time at mid-point, the duration would be:

2(t0-t)/.sol[[1]]

Out[]=

2Log[-1+L]

For the epidemiologic data with t in days, the units are days.

Logistic Distribution

If 1 is substituted for L, the logistic function becomes a distribution function with domain, {-∞, ∞} and range {0, 1}. To put it into the form of the Mathematica LogisticDistribution[], we make some additional substitutions.

-k(t-t0)



/.{L1,t0μ,k1/β,tx}CDF[LogisticDistribution[μ,β],x]

Out[]=

x-μ



Out[]=

x-μ



We have already seen the density function which is the first derivative.

PDF[LogisticDistribution[μ,β],x]

Out[]=

x-μ



1+

x-μ





The differential equation for the distribution function is:

′

(x)

F(x)(1-F(x))

In this form it becomes clear the the term (1 - F(x)) is the survival function of the distribution.

DSolve[{F'[x]F[x](1-F[x])/β,F[μ]1/2},F[x],x]

Out[]=

F[x]

x/β



x/β





This can be simplified.

Exp[-x/β]NumeratorDenominator

x/β



x/β



//Simplify;#1/#2&@@%

Out[]=

-x+μ



A statistical density function is the first derivative of the distribution function. The differential equation in this form is true. The concept of the (1 - F[x]) term as a survival probability is used later to justify the use of the logistic function to model quarantine limited epidemics.

PDF[LogisticDistribution[μ,β],x]==CDF[LogisticDistribution[μ,β],x]SurvivalFunction[LogisticDistribution[μ,β],x]/β//Simplify

Out[]=

True

Other calculations are built in.

Mean@LogisticDistribution[μ,β]Median@LogisticDistribution[μ,β]Variance@LogisticDistribution[μ,β]StandardDeviation@LogisticDistribution[μ,β]

Out[]=

πβ

Hyperbolic functions show up in the characteristic function, which can be obtained eventually by Integrate. Mathematica cannot seem to simplify this without removing the term which adds the location parameter, then putting it back.

cf0=IntegrateExp[Itx]

x-μ



1+

x-μ





,{x,-∞,∞},Assumptions{t∈Reals,x∈Reals,μ∈Reals,β>0}cf0/.μ0//FunctionExpand//Simplify//Exp[Iμt]#&

Out[]=



tβHypergeometric2F11,1-tβ,2-tβ,-





+tβ



tβHypergeometric2F11,1+tβ,2+tβ,-





-tβ

Out[]=

tμ



πtβCsch[πtβ]

The characteristic function is, of course, already built in.

CharacteristicFunction[LogisticDistribution[μ,β],t]

Out[]=

tμ



πtβCsch[πtβ]

When the mean is zero, this is a real function symmetric about the y-axis.

Plot[πtβCsch[πtβ]/.β1,{t,-5,5}]

Out[]=

The density can be recovered easily, from the hyperbolic form.

InverseFourierTransform[

tμ



πtβCsch[πtβ],t,x,FourierParameters{1,1}]PDF[LogisticDistribution[μ,β],x]//ExpToTrig//Simplify

Out[]=

Sech

x-μ

2β



4β

Out[]=

Sech

x-μ

2β



4β

But the forward statistical Fourier transform is not helped by the trigonometric expression.

pdf=FourierTransform

Sech

x-μ

2β



4β

,x,t,FourierParameters{1,1},Assumptions{t∈Reals,x∈Reals,μ∈Reals,β>0};//AbsoluteTimingpdf

Out[]=

{25.4723,Null}

Out[]=

Cosh

-Sinh

Cosh

+

Cosh

-tβHypergeometric2F11,1+tβ,2+tβ,-Cosh

+Sinh

-

Hypergeometric2F11,1+tβ,2+tβ,-Cosh

+Sinh

+Sinh

+

Sinh

-tβ(-+tβ)Hypergeometric2F11,1-tβ,2-tβ,-Cosh

-Sinh

Cosh

2μ

+Sinh

2μ



pdf/.μ0//Simplify//FunctionExpand//ComplexExpand//Simplify//Exp[Iμt]#&

Out[]=

tμ



πtβCsch[πtβ]

The distribution and logistic functions can also be expressed in hyperbolic form.

ExpToTrig@CDF[LogisticDistribution[μ,β],x]ExpToTrig

-k(t-t0)





Out[]=

1+Cosh

-Sinh



Out[]=

1+Cosh[kt-kt0]-Sinh[kt-kt0]

A Closer Look at the Differential Equation

The differential equation for the logistic equation is simplified below.

f'[t]kf[t](1-f[t]/L)

Changing f[t] to x, where x is the number of cumulative cases, the model is shown below.

In[]:=

deModel=kx1-

;

Note that this is a quadratic function. The equation has roots at zero and L, and its peak value occurs when x = L/2.

Solve[deModel0,x]D[deModel,x]//ExpandAllSolve[%0,x]

Out[]=

{{x0},{xL}}

Out[]=

2kx

Out[]=

x



A plot of this quadratic equation is shown below. The epidemic ends when the accumulated case count reaches the limit L. Between 0 and L, the growth rate is always positive and the equation is symmetric about L/2 on the x-axis. It will be shown in the next sections that the epidemiologic data is consistent with this model.

Out[]=

The stream plot shows how the logistic function evolves from the differential equation. The range of the function, below is between 0 and 1000. Along the zero line the vectors diverge and they converge at the 1000 limit.

Out[]=

Why the Model Works

The original derivation of the model was for population growth in a finite space with finite resources for population growth. The idea is that any space can only support some maximum number of inhabitants.

For the model to be applied to an epidemic requires only that there is an accurate and fast method of diagnosis of the infection and that effective quarantine can be applied to the population. Information about the contacts and possibly exposed individuals must also accumulate faster than the infection can spread.

At the onset of an epidemic a small number of cases are identified; an army of public health officials determines the quarantine set, Q, of people who are infected or have been exposed to the people who have the disease. Initially as more infected persons are discovered, Q grows. For an ideal quarantine model within Q people are isolated from the general population and from each other. The goal is that the infection cannot spread beyond Q or within Q, unless there was prior contact with another infected person before the individual was placed in Q. If the quarantine is successful then the set, L, all the people who have been or will eventually become infected, becomes a subset of Q. When that happens no further cases appear outside Q. The exact individuals in Q, who will belong to L, are not knowable initially, but the size of L can in fact be determined rather quickly from the logistic function once L becomes a subset of Q.

The graphic below attempts to illustrate the idea. Click the play button. At time = 0, there is no quarantine, as time increases, the quarantine set, Q, grows faster than the infected set L, once Q contains L, L stops growing, even though its elements may not all yet be fully countable.

Manipulate[ParametricPlot[Evaluate[{Re[#],Im[#]}&/@{If[Time<1/2,(9/10Time+1/10),(9/20+1/10)]Exp[Iθ],11/10TimeExp[Iθ]}],{θ,-Pi,Pi},AxesNone,PlotRange{{-1,1},{-1,1}},PlotLegends{"L","Q"}],{{Time,0},0,3/4,Appearance"Open"},ContinuousActionTrue,SaveDefinitionsTrue]

Out[]=

Now recall the logistic differential equation,

f(t)

t

kf(t)1-

f(t)

If the term,

f(t)

, is removed, the remaining equation is the differential equation for exponential growth. If the quarantine measures don’t succeed and L keeps growing, such that L ≅ r f(t), then the f(t) terms cancel and you are back to exponential growth at a different rate.

f(t)

t

kf(t)1-

f(t)

rf(t)

Eventually the whole population will have the infection, but logistic dynamics may come into play later as the recovered population, which is immune to reinfection, becomes large relative to the remaining susceptible population such that it slows the spread of infection.

With an additional assumption that the probability of the infection arising spontaneously in the population is zero, spread of the infection must be from contact with a person who already is infected, i.e. the infection cannot be transmitted by a non-infected person, then

′

(t)cannotbefasterthankf(t)

and the differential equation,

′

(t)=kf(t)

, gives the upper limit of growth. Note: it is immediatly possible to imagine some rogue health care worker disposing of virus contaminated material in a way that it could infect a huge number of people, so that this assumption could be violated.

In the real world, the cumulative case data at the beginning of the epidemic will be inaccurate. All cases and contacts cannot be discovered instantly; since case data is added when it is discovered, there will be a lag in the cumulative case count as a function of time. If most all of the cases are ultimately discovered, this lag will create a smaller error relative to the current number of cases as time progresses. A consequence of this reality is that data points late in the epidemic should be more accurate. Since it only takes three points to determine a logistic function, a small subset of late points may give a more accurate curve than a data set which includes the counts early in the epidemic. A logarithmic plot of the data is useful to note when growth becomes slower than exponential.

Model Assumptions and Derivation

1. Quarantine measures provide the only method for control of infection spread, such that there can be a maximum number of infected persons, L, which is less than the total susceptible population.

1a. If treatment for the infection or vaccination becomes available during the epidemic, the rate of spread of infection might change as the epidemic progresses. Late in an epidemic that was not controlled, recovered individuals with immunity should also slow the spread. For the onset of a novel viral infection, statement #1 is true, and the other considerations in #1a are false.

2. The maximum growth rate of the number infected is a function of the number infected at that time, i.e.

′

(t)kf(t)

. See the discussion in the section, “Why the Model Works”.

2a. In the real world, this may appear not to be true early in the outbreak, because at the time the outbreak is first noticed, there will be a backlog of infected cases which have not been discovered. They accumulate at the times they are discovered, so it is possible that early data may show exponential or faster growth rates.

3. Information about exposed contacts with infected persons must be collected at faster rates than the infection can spread in the population. Without this assumption, quarantine methods cannot succeed.

4. Quarantine methods can establish a quarantine state for a set, Q, of the population at risk, which completely contains the set, L, of all the people who will be infected by the epidemic outbreak. The set L has a finite size which limits the spread of the epidemic. Refer to section, “Why the Model Works”.

5. The differential equation for exponential growth,

′

(t)kf(t)

, needs a modification for quarantine methods. In a quarantine an instantaneous transmission of the virus from infected person, A, to a non-infected person B, may not happen. So the equation is modified.

′

(t)kf(t)P(t)

Where P(t) is an instantaneous probability at time, t, that the transmission will occur. The range of P(t) is between 0 and 1. Given the assumptions above, P(t) is the probability generated by a survival function,

P(t)1-

f(t)

where L exists as a limit. This multiplier, P(t), functionally reduces k, the growth rate, toward zero as f(t) approaches L
Completing the substitution yields the logistic differential equation.

′

(t)kf(t)1-

f(t)

6. In the real world the the exponential growth of the infected population and the quarantine operations are going on simultaneously and competing with each other. Initially the quarantine is incomplete. As cases are put into the quarantine set, Q, the growth rate of the infection may be slowed, but it will not have a limit. Once the quarantine set, Q, fully contains L, all the people who can be infected, the epidemic theoretically is guaranteed to end in the future, and new cases should only arise within Q. Cases in the set, L, become countable only when a diagnosis is made by laboratory test or clinical conditions. The diagnosis cannot be made until after an incubation period, after which symptoms develop and virus particles are released from cells into secretions.

The accumulated case count, f(t), given in the data from an epidemic represents the information flow from the epidemic at time, t. Thus f(t) lags the actual state of the epidemic at time, t, because of the delay in diagnosis. Once the set L is fully contained by Q, the growth rate is fully controlled by the logistic function rather than the exponential function. This is true even if L cannot yet be fully counted. The case counts are cumulative; after L is contained, the data points should be expected to more accurately fit the logistic function.

Application of the Logistic Model to Actual Data

Procedures

In[]:=

Clear[modelFitDisplay];modelFitDisplay[NLM_,data_,opts___]:=Module[{bands,p0,p1},bands[t_]=NLM["SinglePredictionBands",ConfidenceLevel0.95];p0=Plot[{NLM[t],bands[t]},{t,1,data[[-1]][[1]]+5},opts,Filling{2{1}},PlotStyle{Darker@Red,Darker@Orange,Darker@Orange},GridLinesAutomatic,PlotRangeAll,opts];p1=ListPlot[data,PlotRangeAll];Show[p0,p1]];

In[]:=

Clear[derivativeDisplay];derivativeDisplay[NLM_,data_,opts___]:=Module[{daynum,end,peak},daynum=data[[-1]][[1]];peak=t0/.NLM["BestFitParameters"];end=x/.FindRoot[NLM'[x]1,{x,peak,peak+200},(*forcetorighttail*)Method"Brent"][[1]];Print["Last case: ",DateString[DatePlus[{2020,1,21},Ceiling@end],"DateShort"]];Show[Plot[{NLM'[x],NLM''[x]},{x,data[[1]][[1]]-4,daynum+5},opts,GridLinesAutomatic,PlotRangeAll,PlotLegends{"First Derivative","Second Derivative"}],ListPlot[Tooltip[#,#]&/@{{daynum,NLM'[daynum]}},PlotStyle{Red,PointSize[0.02]}]]];

In[]:=

Clear[dataD1];dataD1[data_,{k_,L_}]:=Table[{data[[i]][[1]],kdata[[i]][[2]](1-data[[i]][[2]]/L)},{i,Length[data]}];

In[]:=

logistic=L/(1+Exp[-k(t-t0)]);

In[]:=

Clear[logisticDEFit];logisticDEFit[data_,k0_:0.3,L0_:2000,graphic_:True,optsInterpolation___]:=Module{dataFn,dataDE,dataDEPlot,deModel,deNLM,bands},dataFn=Interpolation[data,optsInterpolation,Method"Spline"];dataDE={dataFn[#],dataFn'[#]}&/@data[[All,1]];dataDEPlot=ListPlot[dataDE,PlotStyleDarker@Red,GridLinesAutomatic,PlotRangeAll];deModel=kx1-

;deNLM=NonlinearModelFit[dataDE,{deModel(*,L>0*)},{{k,k0},{L,L0}},x];If[graphic,bands[L_]=deNLM["SinglePredictionBands",ConfidenceLevel0.95];Print@Show[Plot[{deNLM[x],bands[x]},{x,0,Abs@L/.deNLM@"BestFitParameters"},GridLinesAutomatic,Filling{2->{1}},FrameTrue,FrameLabel{"Cases","Cases Per Day"},PlotRangeAll,PlotLabel"Logistic Differential Equation"],dataDEPlot]];Flatten@{deNLM@"BestFitParameters",FindRoot[dataFn[t0]L/2/.deNLM["BestFitParameters"],{t0,data[[-1]][[1]]}]};

In[]:=

Clear[logisticParmProg];logisticParmProg[data_]:=Module{nlmd,data1},nlmd=NonlinearModelFit#,

-k(t-t0)



,{{k,0.3},{L,5000},{t0,50}},t&/@Table[Take[data,i],{i,4,Length@data}];data1=Drop[data[[All,1]],3];GraphicsColumn[{ListPlot[Transpose[{data1,L/.#["BestFitParameters"]&/@nlmd}],PlotRangeAll,GridLinesAutomatic,FrameTrue,FrameLabel{None,"L Parameter"}],ListPlot[Transpose[{data1,k/.#["BestFitParameters"]&/@nlmd}],GridLinesAutomatic,FrameTrue,FrameLabel{None,"k Parameter"},PlotRangeAll],ListPlot[Transpose[{data1,t0/.#["BestFitParameters"]&/@nlmd}],GridLinesAutomatic,FrameTrue,FrameLabel{"Day Number","t0 Parameter"},PlotRangeAll]},Spacings0];

In[]:=

epidemicDates[NLM_FittedModel]:=DatePlus{2020,1,21},Floort0-

Log[-1+L]

/.NLM@"BestFitParameters",DatePlus{2020,1,21},Ceilingt0+

Log[-1+L]

/.NLM@"BestFitParameters";

The data below are derived from the JHU CSSE database for South Korea. The day numbers on the abscissa are based on day one being 22 January 2020, which is the first data day of database. The first 29 days from the data represented 31 accumulated cases, which had not shown any growth, these were dropped and 31 was subtracted from the accumulated cases. The data thus represent the outbreak which began in South Korea on February 19, 2020, and showed up in the database the next day. The logarithmic plot below shows the data becoming concave downward after the first three points, which show exponential growth.

Out[]=

The NonlinearModelFit to the logistic function is shown below. The graph is shown with 95% confidence interval along with the data points.

skoreaNLM1=NonlinearModelFit[sk2,logistic,{{k,0.3},{L,10000},{t0,50}},t];skparm=skoreaNLM1["BestFitParameters"]modelFitDisplay[skoreaNLM1,sk2,FrameTrue,FrameLabel{"Day Number","Cumulative Cases"},PlotLabel"South Korea Cases",PlotRange{{sk2[[1]][[1]]-5,sk2[[-1]][[1]]+5},All}]

The next graph shows the first and second derivatives of the logistic function above. The peak of the epidemic occurs when the second derivative curve crosses zero. The green dots represent points calculated from the data points, using the differential equation,

′

(t)kf(t)1-

f(t)

, and the values for k and L determined by the fit. The red dot is the location on the first derivative curve at the last day of the data. The last case date is the day when the value of the first derivative becomes less than one.

kdp1=derivativeDisplay[skoreaNLM1,sk2,FrameTrue,FrameLabel{"Day Number","Cases per Day"},PlotLabel"South Korea"];Show[kdp1,ListPlot[dataD1[sk2,{k,L}/.skoreaNLM1["BestFitParameters"]],PlotStyleDarker@Green]]

The graphic below shows parameters calculated by successive parameter fits, daily, starting with day 4 of the data set. By day number 43, after just 14 data points, the parameters begin to show convergence to a stable solution. At this point, L can be reasonably estimated and the epidemic should end as the number of cases approaches L. The convergence, to a stable logistic solution suggests the the quarantine measures have succeeded and the epidemic is in the logistic phase in South Korea.

Out[]=

With the parameters converging, the duration of the epidemic can be calculated in days.

2Log[-1+L]

/.skoreaNLM1@"BestFitParameters"//Ceiling

Out[]=

Estimated epidemic start and end dates from the logistic function equation, rather than the derivative equation, which was used for last case date.

skDates=epidemicDates[skoreaNLM1]

Out[]=

{{2020,2,4},{2020,3,29}}

Here is another potentially useful graphic. This shows the daily differences of the natural logarithm of the daily case number. The actual log differences from the data are in orange and the predicted values from the parameter fit to the logistic model are in blue. With the epidemic in control they are clearly converging. With exponential growth, the orange points would be in a horizontal line at the growth rate. Bad data points tend to vary widely.

Out[]=

	Predicted
	Actual

The data can also be fit to the model using the logistic differential equation. To do this a method is needed to calculate the first derivative from the data. In the graphics below the Interpolation function is used with InterpolationOrder  {1, 2, 3} sequentially. The default Method  “Spline” is used to to give a continuous derivative function. The red dots are the first derivative points calculated at the data point; note that they change with different interpolation orders so that the method may not be very precise.

Print[logisticDEFit[sk2,0.3,9000,True,InterpolationOrder#]]&/@Range[3];

The South Korean data are rather ideal compared to most of the data sets. The log plot of the early data from Washington state is shown below. The first five points seem to be growing exponentially; the next five points also seem exponential at a slower rate. Finally the last five points begin to show a pattern consistent with the logistic model.

Out[]=

This could be a case of data collection lagging at the onset of the outbreak, with old cases being discovered at an exponential rate, well after they actually occurred. The next graph shows the fit to the logistic model using all 15 data points.

waNLM=NonlinearModelFit[wa,logistic,{{k,0.3},{L,1000},{t0,60}},t];waNLM@"BestFitParameters"modelFitDisplay[waNLM,wa,PlotRange{{wa[[1]][[1]]-2,wa[[-1]][[1]]+5},All},PlotLabel"Washington Fit 15 Points",FrameTrue,FrameLabel{"Day Number","Cases"}]

When the fit is made with the last five points, the L parameter is lower. It is reasonable to expect more accurate parameter estimates going forward, if analysis starts with the later points, which contain all the cumulative information since the onset. At this point it is too soon to make reasonable projections.

waNLM=NonlinearModelFit[Take[wa,-5],logistic,{{k,0.3},{L,1000},{t0,60}},t];waNLM@"BestFitParameters"modelFitDisplay[waNLM,Take[wa,-5],PlotRange{{wa[[1]][[1]]-2,wa[[-1]][[1]]+5},All},PlotLabel"Washington Fit Last 5 Points",FrameTrue,FrameLabel{"Day Number","Cases"}]

The log plot of the Italian outbreak is below.

Out[]=

The differences of daily Log[cases], model and actual data are below. At this point Italy may be just starting to control the epidemic, a few more data points are needed to establish this.

itNLM=NonlinearModelFit[it,logistic,{{k,0.2},{L,25000},{t0,70}},t];itDates=epidemicDates[itNLM]dtlist=DateRange[DatePlus[{2020,1,21},it[[2]][[1]]],itDates[[2]]];DateListPlot[{Transpose[{dtlist,Differences@Log[itNLM[#]&/@Range[it[[1]][[1]],DateDifference[{2020,1,21},itDates[[2]]][[1]]]]}],Transpose[{Take[dtlist,Length@it-1],Differences@Log@it[[All,2]]}]},JoinedFalse,FrameTrue,FrameLabel{None,"Daily Differences Log[Cases]"},PlotLegends{"Predicted","Actual"},GridLinesAutomatic,PlotLabel"Italy"]

The fit shows a very wide range for the estimate of the parameter, L, the ending number of infected cases.

itNLM@"ParameterConfidenceIntervalTable"

Out[]=

	Estimate	Standard Error	Confidence Interval
k	0.218849	0.0114052	{0.195255,0.242442}
L	76702.8	13294.4	{49201.2,104204.}
t0	57.5663	1.35814	{54.7567,60.3758}

POSTED BY: Robert Rimmer

6 Replies

Sort By: