Introduction

In the section on linear regression, we saw that it is relatively simple to fit a set of data of the form ((Subscript[x, i]),(Subscript[y, i])), i = 1, 2, ..., n to an approximate function of the form

g(x)=(∑_(j=1)^m)(Subscript[a, j])(Subscript[Φ, j])(x).

Requiring that the sum of squared errors be minimal will always give a linear system of m equations for the "weights" (Subscript[a, j]) that can be solved using standard linear algebra techniques.  Since the algorithm is so striaghforward, it has already been implemented in Mathematica, in the Fit command.

A big part of doing linear regression is coming up with a more or less inspired guess for the "basis" functions, (Subscript[Φ, j])(x).  However, it is often the case that the functions we might pick have implicit parameters.  For example,consider the following set of data:

datapoints=Import["T:\\white\\classes\\Computer Aided Math\\ClassNotes\\Mathematica 7\\nlreg1.txt","Table"];

dim=Dimensions[datapoints]

datapointsplot=ListPlot[datapoints,PlotStyle->Black,PlotMarkers->{"o",Small}]

From the graph we can see that the data goes from x = 0 to x = 10, is sinusoidal and that it has a decreasing amplitude.  We can also see that the data does not form a function (look at x = 0 for instance).  This might be two repetitions of an experiment that has a slightly random component.  In order to describe this data with a function we might guess that the periodic component looks like  for some value of k.  To explain the decrease in amplitude we might think of exponential decay.  Hence, a reasonable model for this data could be

g(x)=a (e^(λ x)) sin(k x)

The trouble with this model is that we have no idea what to choose for the parameters k and λ, although if we did know these values we could use linear regression to find a.

The goal of this chapter will be to figure out some ways of estimating parameters which are imbedded in the approximating function, g(x).  One technique is to use nonlinear regression which will be discussed later.  Other methods depend on knowing more knowledge about the data, for instance, that it should have been produced by a differential equation, or have some particular functional form.  These lead to other ways to estimate parameters, which is what we will discuss first.  The data we will use (not the above data) comes from the 1906 plague epidemic in Bombay, studied by Kermack and McKendrick, and the model we will attempt to fit comes from the theory of epidemics.

Basic Theory of Epidemics

In this section we will discuss an example from the theory of epidemics.  In the theory of epidemics, populations of people are divided into three classes:

1.  The susceptible class, S, who have not yet been infected,

2.  The infected class, I, who have caught the disease, and

3.  The  removed class, R, who have recovered (if that is possible) or died.

If N is the total population, then the following relation holds:

N = S + I + R .

This is a statement that these three classes make up the whole population, and there are no others.  For epidemics that occur rapidly, one may neglect the effects of birth rate, natural death, and population redistribution to derive a system of equations relating the three types of people in the population,

(dS/dt)=-β S I ,

(dI/dt)=β S I-ν I ,

(dR/dt)=ν I.

The parameter β defineds a frequency at which encounters between infected individuals and susceptible individuals results in infection, and ν defines the rate at which infection results in death (or recovery).  The reason that the varying terms on the right hand side of the equations are repeated, but with different signs, is that a loss from one equation must be a gain for another.  If a susceptible person becomes infected, then the infected population must grow, and the same with the death of infected people.  This fact is reflected mathematically in the following way.  If you sum all three equations, the right hand sides cancel, giving

(d/dt)(S+I+R)=0.

Integrating this simple equation gives that the total population is constant, as stated above.

The total population, N, is called a "constant of the motion," and can help us solve the equations governing an epidemic.  By replacing the variable S according to the rule

S = N - I - R ,

we can reduce the equations to a system of two equations,

(dR/dt)=ν I.

(dI/dt)=β  I (N-I-R)-ν I ,

Finally, we can differentiate the first equation and substitute into the second:

(((d^2)R)/(dt^2))=ν (dI/dt), and I=(1/ν) (dR/dt),

so we get,

(((d^2)R)/(dt^2))=(β N-ν) (dR/dt)-β R (dR/dt)-(β/ν)((((dR/dt)))^2).

An interesting note in the theory of immunology: if

β N < ν ,

then the right hand side is negative-definite whenever the death rate goes up.  Hence, the number of deaths must decelerate and the epidemic can not become self-sustaining.

We can take a page out of physics and use an energy method to rewrite the equation as

-(β/ν)((((dR/dt)))^2)=(d/dt)((dR/dt)+(1/2)β (R^2)-(β N-ν) R).

If we make the assumption that the left hand side (in particular the derivative of R with respect to t) is small in relation to the other terms (which an immunologist would say means that the rate of removal is small in relation to the removed population), then the above equation can be integrated as follows:

(dR/dt)+(1/2)β (R^2)-(β N-ν) R=K, for some constant K.

Before the epidemic begins, both the rate of removal and the removed population must be zero, so the constant of integration, K, must be zero.  Hence, under simplifying assumptions, the equation for the removed population becomes

(dR/dt)+(1/2)β (R^2)-(β N-ν) R=0 .

This looks very like the logistic equation we met in the sections for population models.  We can put the equation into Mathematica and solve it:

eqn1=DSolve[R'[t]+β R[t]^2/2-(β n-ν) R[t]==0,R[t],t]

If we change variables, letting λ = β N - ν, this solution can be written as

R=((2λ)/(β+(Subscript[C, 1]) λ Exp[-2λ t])).

A different way to write this solution would be

R=(λ/β)(1+Tanh[λ(t-(Subscript[t, 0]))]).

The rate of removal, or number of deaths per unit time would be

(dR/dt)=((λ^2)/β)((Sech[λ(t-(Subscript[t, 0]))])^2).

Curve Fitting by Inspiration

In this section we will examine how we might use our knowledge about the theory of epidemics to fit a functional form to the Bombay plague data.  We will not worry about anything so sophisticated as minimizing error; instead we will just attempt a reasonable guess for parameter values.

The city of Bombay, India kept records of mortality through the plague of 1906.  Over a period of thirty weeks many thousand people died, and on a weakly basis the number of deaths were recorded.  The data is as follows:

dpw={20.,22.,25.,31.,40.,60.,66.,93.,120.,182.,293.,388.,451.,636.,765.,778.,721.,701.,872.,920.,794.,582.,411.,207.,102.,66.,49.,45.,36.,30.};
len=Length[dpw]

wks=Table[j,{j,1,len}]

xypairs=Table[{wks[[j]],dpw[[j]]},{j,1,len}]

dpts=ListPlot[xypairs,PlotStyle->Black,PlotMarkers->{"o",Small},PlotLabel->"Deaths/week"]

Our goal is to fit this data with a functional form and thereby estimate parameters in the differential equations that describe how epidemics proceed.  Notice that the information plotted is in deaths per week, which we would think of as the derivative of R, or the rate of removal.  Hence, we will use inspirational guessing techniques to fit the solution

(dR/dt)=((λ^2)/β)((Sech[λ(t-(Subscript[t, 0]))])^2).

to the plague data.

First notice that the rate of removal is symmetric about the point in time t=(Subscript[t, 0]).  This point is the center of mass of the analytic description for the rate of removal.  It seems reasonable that it should be chosen to match the center of mass of the data.  The center of mass of the data is

t0=Sum[j*dpw[[j]],{j,1,len}]/Sum[dpw[[j]],{j,1,len}]

So, 17.2557 weeks is the point in time about which the epidemic is "centered".

If we had a reasonable value for λ, we could perform a linear regression against the function

(sech^2)(λ(t-(Subscript[t, 0])))

to fit the data.  Recall that the hyperbolic secant is defined in terms of the natural exponential function as follows:

Sech[λ(t-(Subscript[t, 0]))]=(2/(Exp[λ(t-(Subscript[t, 0]))]+Exp[-λ(t-(Subscript[t, 0]))])).

Taking λ > 0, as t-(Subscript[t, 0])→∞, (sech^2)(λ(t-(Subscript[t, 0])))→4Exp[-2λ(t-(Subscript[t, 0]))] .  Thus, if we fit the tail end of the plague data to an exponential function, we may be able to extimate λ.  To fit something exponentially, we wish to find coefficients a and b so that

Exp[a (Subscript[t, j])+b]=(Subscript[((dR/dt)), (Subscript[t, j])]),

where we may take the equals to mean a least-squares best fit and the right hand side is the experimentally measured removal rate at week (Subscript[t, j]).  Hence, taking the logarithm of the data, we would have

a (Subscript[t, j])+b=ln((Subscript[((dR/dt)), (Subscript[t, j])])),

meaning that we wish to find a line which fits the logarithm of the data according to least squares.  This isn't too hard in Mathematica.  The commands are:

tail1=dpw[[1;;8]]
tail2=dpw[[23;;30]]

Take the first and last eight data points on the removal rate, that are the visual "tails" of the data and isolate them from the bulk of the data.  We can put the logs of the tails together:

tail=Join[Log[tail1],Log[tail2]]

the t values that go along with this list is:

twks=Join[wks[[1;;8]],31-wks[[23;;30]]]

then fit the data

logfit=Fit[Table[{twks[[j]],tail[[j]]},{j,1,16}],{1,x},x]

The coefficient of x in this expression is twice the λ we are looking for, so

λ=D[logfit,x]/2

Now we can do a linear regression

eq=Fit[xypairs,{Sech[λ*(t-t0)]^2},t]

Graphing this equation along with the data gives:

Show[
dpts,
Plot[eq,{t,0,30},PlotStyle->Red]
]

This fitted function seems to do a fair job of capturing the initial growth of the data, and peaks out more or less where the actual plague data peaks.  It misses the complicated behavior at the peak of the death rate distribution, and doesn't seem to do very well as the epidemic dies out.  On the other hand, it's not clear that the weird behavior at the peak of the epidemic is not due to some effect we are not accounting for.  For example, people trying hard not to come in contact with plague victims at the peak of the epidemic, or perhaps a two-week spell of cooler weather.  The original researchers, following pretty much the methodology we have, predicted an amplitude of 890, λ = 0.2, and (Subscript[t, 0])=17.  Given the fact that both sets of researchers, we and they, neglected a term in our simplified model, maybe we've done as well as can be expected with this methodology.

Nonlinear Regression

An alternative to our heuristic fitting technique would be to attempt a full nonlinear regression.  As with linear regression, the idea is to minimize the square error between the actual data and some representation we have chosen.  Let g(t,(Subscript[a, 1]),(Subscript[a, 2]),...,(Subscript[a, m])) be a model representation we have chosen for the data, that depends on the parameters (Subscript[a, 1]), (Subscript[a, 2]), ..., (Subscript[a, m]).  If the data are given by the set of points ((Subscript[t, j]), (Subscript[f, j])), for j = 1, ..., n, then the error can be defined as

E=(∑_(j=1)^n)((((Subscript[f, j])-g((Subscript[t, j]))))^2),

where E is a function of the parameters (Subscript[a, k]), for k = 1, 2, ..., m.  As before, our goal is to minimize the size of the error.  From multivariable calculus, we know that at such a minimum the partial derivatives must be zero, or

((∂E)/(∂(Subscript[a, k])))=-2((∑_(j=1)^n)((∂/(∂(Subscript[a, k])))g((Subscript[t, j])))((Subscript[f, j])-g((Subscript[t, j]))))=0.

This gives a system of m equations for the parameters, (Subscript[a, k]), at which such a minimum can occur.  The problem is that the equations are essentially nonlinear.  In general, there is no way to solve the nonlinear equations.  However, since we are already dealing with numerical data, looking for a numerical approximation to the solution is sufficient.

Let's see how this works with our plague problem.  In this case the representative function has the form:

g(t,a,λ,(Subscript[t, 0]))=a (sech^2)(λ(t-(Subscript[t, 0]))),

and the error function is deined above.  In Mathematica we can accomplish this as follows.  The function we are fitting is

Clear[λ,t0,a]

g[t_]=a*Sech[λ*(t-t0)]^2

The error function would be

err=Sum[(dpw[[j]]-g[wks[[j]]])^2,{j,1,len}];

The three partial derivatives of this error function are given by

eqna=D[err,a];
eqnl=D[err,λ];
eqnt0=D[err,t0];

Next we want to solve for a critical point.  We find a numerical approximation for the critical point of interest by using the NSolve command.  Since we are trying to solve complicated nonlinear equations the numerical routine is much more efficient if we can give it a good initial guess for the solution.  Based on the work we did above, a good initial guess might be a = 890, λ = 0.2, and t0 = 17.  Plugging these values into the NSolve command gives:

cpts=FindRoot[{eqna==0,eqnl==0,eqnt0==0},{{a,890},{λ,0.2},{t0,17}}]

To plot this new approximating function we would do the following:

nlp=Plot[g[t]/.cpts,{t,0,30},PlotStyle->Blue]

Combining this with the data point gives

Show[dpts,nlp]

Compare this with the result we obtained by intuition:

Show[
dpts,nlp,
Plot[eq,{t,0,30},PlotStyle->Red]
]

It is clear that the nonlinear regression approach does a much better job of approximating the data then the educated guessing we did earlier.  However, it still misses the complicated behavior near the peak and does not do a great job at the end of the epidemic.  As before, given that we have neglected terms in an attempt to simplify the differential equation used to model the epidemic, we've done reasonably well.

Fitting the Differential Equation

The trouble with linear regression is that it can only give the relative weights of a linear combination of known function in a fit.  Nonlinear regression is more general, in that it can be use to find parameters wherever they occur in an approximating function.  However, in each case you have to have a guess as to the form of the approximating function or functions.  In the epidemic case we were provided with a reasonable guess for the basis function: a square hyperbolic secant function.  We got this guess by making some simplifying assumptions about the differential equation that models the epidemic so that we could solve the DE in closed form.  Instead of using this approximation, we might consider fitting the original DE to the data and then solve the DE numerically.

Recall the equation of interest for the number of deaths, R,

(((d^2)R)/(dt^2))=(β N-ν) (dR/dt)-β R (dR/dt)-(β/ν)((((dR/dt)))^2).

If we again write λ = β N - ν, then we get the equation:

(((d^2)R)/(dt^2))=λ (dR/dt)-β R (dR/dt)-(β/ν)((((dR/dt)))^2).

If we knew λ, β, andν, then we could solve this equation numerically and fit that solution to the data boundary conditions.  The right hand side of the above expression is a polynomial in R and (dR/dt), and we know from experience that polynomials are pretty easy to fit.

To make things a little more clear, let's change variables by letting

x = R, y=(dR/dt), and f(x,y)=(((d^2)R)/(dt^2)).

Then the differential equation becomes a function of two variables,

f(x,y)=λ y-β x y-(β/ν)(y^2).

The data we have on the plague in Bombay then needs to be in terms of ((Subscript[x, j]),(Subscript[y, j])), which is the total number of deaths at week j and the deaths per week at week j, respectively.  We do not have direct data on the total number of deaths, (Subscript[x, j]), nor do we have direct data for f((Subscript[x, j]),(Subscript[y, j])), the derivative of the removal rate at week j.  However, we generate this data based on what we do have.

The total number of deaths is just the sum of the weekly deaths over the first j weeks.  That is,

(Subscript[x, j])=(∑_(k=1)^j)(Subscript[y, k]).

The dependent variable, f, is really the derivative in time of the rate of removal, y.  From basic calculus, we know that (using Mathematica's notation)

(du/dt)=Limit[((y(t+h)-y(t-h))/(2h)), h=0],

provided the limit exits.  Thus, we may make the following approximation for (Subscript[f, j])=f((Subscript[x, j]), (Subscript[y, j])):

(Subscript[f, j])≈(((Subscript[y, j+1])-(Subscript[y, j-1]))/2).

Note that the denominator stands for two weeks.  We should note that this approximation is only valid for j = 2, 3, ..., n-1.  In Mathematica, we can generate (Subscript[x, j]), (Subscript[y, j]), and (Subscript[f, j])for j = 2, 3, ..., n-1, as follows

Clear[x,y,f,λ,β,ν]
DEdata=Table[{Sum[dpw[[k]],{k,1,j}],dpw[[j]],(dpw[[j+1]]-dpw[[j-1]])/2},{j,2,len-1}];

Now we want to do a best (linear) fit to the polynomial above:

eqde=Fit[DEdata,{y,x*y,y^2},{x,y}]

Then the parameters we are after would be

λ=Coefficient[eqde,y]/.{x->0}

β=-Coefficient[eqde,x*y]

ν=-β/Coefficient[eqde,y^2]

Now we can use these values to solve the differential equation numerically.  Note that we will need initial conditions to solve the DE.  We can use any one week as our initial conditions.  Also note that it may be better to convert to a system of first order differential equations:

(dR/dt)=U,

(dU/dt)=λ U-β R U-(β/ν)(U^2).

To solve these equations we would do the following:

de={r'[t]==u[t],u'[t]==λ*u[t]-β*r[t]*u[t]-β/ν*(u[t])^2}

sols=NDSolve[Join[de,{r[20]==DEdata[[19,1]],u[20]==DEdata[[19,2]]}],{r,u},{t,0,30}]

ssp=Plot[Evaluate[u[t]/.sols],{t,0,30},PlotStyle->Darker[Green]];
Show[ssp,dpts]

Notice that we have used week 20 as our initial condition.  If we want to use week 5 instead, we would get:

sols=NDSolve[Join[de,{r[5]==DEdata[[4,1]],u[5]==DEdata[[4,2]]}],{r,u},{t,0,30}];
ssp=Plot[Evaluate[u[t]/.sols],{t,0,30},PlotStyle->Darker[Green]];
Show[ssp,dpts,PlotRange->All]

Clearly the choice of initial condition will effect the outcome of this method.  We can now experiment with the fits generated by using various initial conditions from the data.  The fit generated by the week 20 data may not be more accurate then the nonlinear regression we did before, but it does have a more qualitative correctness.  That is, including all of the terms in the differential equations captures trends in the data more accurately.  For example, the fact that the epidemic dies out more rapidly than it starts up.

None of the techniques for parametric estimation capture the strange behavior of the data near the peak of the epidemic, which is part and parcel for trying to model "real life" data with relatively simple mathematical models.

Exercise 9

In earlier sections, we examined solutions to some different population-dynamic models.  The most realistic model for a single population with limiting growth due to resource consumption was the logistic model,

(dP/dt)=a P (1-b P).

A great deal of data has been taken to see whether or not this theoretical model is valid.  Bodenheimer and Nolan took the following growth data on bee hives near Baltimore (there should be 19 data points):

Days={7,14,21,28,35,42,49,56,63,70,77,84,91,98,105,112,119,126,133};
Bees={0.6,1.0,1.4,2.0,2.5,3.3,6.8,10.5,15.9,20.2,30.0,39.7,47.2,54.6,61.4,65.9,71.4,76.4,76.9};

where the number of bees is counted in thousands.  Below you will explore how well the logistic equation describes this data set.  Be clear about what assumptions you make and your rationale for making them.

1.

Approximate the derivatives of the population and choose parameters to fit the logistic equation to the data.  How well do solutions to the fitted equation fit the data?

2.

Choose parameters in the logistic equation to fit the asymptotic behavior of the data.  Fit an analytic solution with these parameters to the data by choosing the free constant to match something else you know about the behavior of the logistic equation.  Compare with the results from problem 1.

3.

Use nonlinear regression to fit logistic equation solutions to the data.  Compare your fit with the data and with the results from problems 1 and 2.

4.

Based on your work with this data set, is the logistic equation a good description of the growth of a bee hive?  What are the pros and cons of the three methodologies applied above?  In what circumstances is it appropriate or inappropriate to apply each methodology?