Curve fitting--4

S. Arlinghaus
Logistic curve fitting.
 
     The exponential function leads to a situation of unbounded population
growth.

Assumption for the exponential:
     The rate of population growth or decay at any time t is proportional
to the
size of the population at t.

Symbolically:  let Y_t  [N from last time] represent the size of a
population at time t.
The rate of growth of Y_t is proportional to Y_t:
        
dY_t / dt = kY_t

where k is a constant of proportionality.

This assumption led to the curve Y_t = Y_(t0) e^(kt) ---- 
a situation in which growth is unbounded as t becomes large.

Bounded growth--the logistic function.

     Assume now that in reality, when the population gets large,
environmental factors
dampen growth ---  a population-environment dynamic.

To express this idea symbolically:
     assume that eventually the growth rate decreases ---- dY_t / dt
eventually decreases.
     assume that population size is limited to some maximum, q, 
          where 0 < Y_t < q.  As Y_t approaches q, dY_t / dt approaches 0,
so that population
          size tends to become stable as t approaches infinity.
The result is a model that is exponential in shape initially and includes
effects of
environmental resistance in larger populations.  A typical geometric
characterization is             
shown below.            



                










                                        t

The various coordinates in this figure are derived from expressing the
assumptions of bounded                                  
growth algebraically.  One such expression is:                                  

dY_t / dt = kY_t * (q-Y_t)/q                                    

This equation meets the conditions of bounded growth because in the factor
(q-Y_t)/q                                       
when Y_t is small                                       
(q-Y_t)/q is close to 1 and the growth is therefore close to the
exponential in the early stages;                                        
when Y_t is large--that is, close to q--                                        
(q-Y_t)/q is close to 0, and the growth rate dy_t / dt tapers off to the
horizontal asymptote.                                   
The factor (q-Y_t)/q, which is what alters the symbolic form of the
logistic from the exponential,                                  
acts as a damper to growth--as required.                                        
*************************************************************************************************************                                   
To understand where the various quatities in the graph above come from,
consider the following                                  
set of manipulations.

Replace k/q by K so that the logistic equation now reads:

dY_t / dt = KY_t(q-Y_t)

in which the rate of growth is proportional to the product of the
population size and the
 difference between the maximum size and the population size.

Solve this latter equation for Y_t  ----  separate the variables of the
differential equation

(dY_t)/(Y_t(q-Y_t)) = K dt

to yield

INT (dY_t)/(Y_t(q-Y_t)) = INT K dt.

Use a Table of INtegrals on the rational form in the left-hand integral:

1/q * ln|(Y_t)/(q-Y_t)| = Kt +C or,
ln|Y_t/(q-Y_t)| = qKt +qC.

Because Y_t > 0 and q-Y_t > 0,

ln (Y_t/(q-Y_t)) = qKt +qC.

Therefore,

Y_t/(q-Y_t) = e^(qKt)e^(qC).

Replace e^qC by A.  Therefore

Y_t/(q-Y_t) = Ae^(qKt)
Y_t=(q-Y_t)Ae^(qKt)
Y_t=qAe^(qKt)-Y_tAe^(qKt)
Y_t+Y_tAe(qKt)=qAe^(qKt)
Y_t = (qAe^(qKt)/(Ae^(qKt) + 1))

Now divide top and bottom of the last equation by Ae^(qKt) [equivalent to
multiplying by 1],
so that

Y_t = q/(1+(1/(Ae^(qKt)))) = q/(1+ 1/A * e^(-qKt)).

Replace 1/A by a and -qK by b producing a common form for the logistic
function:

Y_t = q/ (1 + ae^(bt))

with b < 0 because b = -qK, and q, K > 0.
*********************************************************************************************************

Facts about the graph of the logistic equation.

a.  The line Y_t = q is a horizontal asymptote for the graph.

     This is so because, for b<0,

lim (t --> inf) q/(1+ae^(bt)) ------> q/(1+a(0)) = q.

Can the logistic curve cross this horizontal asymptote to produce one or
multiple crossings?

















That is, can the curve approach the asymptote from above, or might the
curve oscillate back
and forth across the asymptote while closing down on it?

This idea can be expressed in notation by setting Y_t = q in the logistic
equation and solving
to see if there are any points that lie simultaneously on the curve and
the asymptote.

Can it be that:

Y_t = Y_t / (1+ae^(bt))?

Or,

1=1+ae^(bt)?

Or, ae^(bt) =0--or that a=0--no, becasue a=1/A.  Or that e^(bt)=0--no.

Thus, the logistic growth curve cannot cross the horizontal asymptote so
that it approaches
it entirely from one side--in this case, from below.

b.  Find the coordinates of the inflection point of the logistic curve.

The inflection point is the point at which the curve changes from concave
up to concave down.
It is also the point at which the growth rate is a maximum--to find this
maximum differentiate one
form of the logistic equation:

logistic equation:  dY_t /dt = KY_t(q-Y_t)=KqY_t - K(Y_t)^2

derivative:  d^2Y_t /dt^2 = Kq - 2KY_t

Set this last equation equal to zero.

Kq - 2KY_t = 0.

Therefore, Y_t = q/2.  This is the vertical coordinate of the inflection
point of the curve for
Y_t, the logistic curve  --  dY_t /dt is increasing to the left of q/2 and
increasing to the right
(from evidence of second derivatives).   The value found is a maximum.

Horizontal component:
     To find t, put Y_t = q/2 in the logistic equation and solve:

q/2 = q/(1 + ae^(bt))

Solving, 1 + ae^(bt) = 2; e^(bt)=1/a; e^(-bt)=a; -bt=ln a, and

t = (ln a) / (-b)

Thus, the labelling of the logistic curve is as given earlier--
the height of the inflection point is half the height of the horizontal
asymptote.

LOGISTIC CURVE
FIT USING ACTUAL DATA.

EQUATION:
     y=q/(1+ae^(bt)), b<0.
 
FIND:
     from the given data, values for a and b.
     Then, use the equation on a number of  t-values to graph
     the equation.

TO FIND a and b:

     All that is required is two values of  t, an upper bound q,
     and the assumption that the distribution is to be logistic.

EXAMPLE:
     Inland coal consumption in India--'000 tons
Data of K. Hornbarger (545 student in 1992).

t       y       least sqs
year 0 is               exponential fit
1970; year 16           
is 1987         
projected beyond 1987           
0       77435   74332.97938 
1       81483   78484.54458 
2       82077   82867.97851 
3       87569   87496.23125 
4       92759   92382.97616 
5       102426  97542.65027 
6       105436  102990.4969 
7       101409  108742.6108 
8       106165  114815.9855 
9       108040  121228.5638 
10      127545  127999.2905 
11      136160  135148.1685 
12      145569  142696.3179 
13      148470  150666.0384 
14      163950  159080.8751 
15      174324  167965.6881 
16      184790  177346.7261 


Now suppose instead that estimates based on remaining reserves          
require that inland consumption never exceed 250000 tons of coal                
per year.  Suppose also that the last year for which we have actual             
data is the year 1987.  Assume a logistic representation in             
the form y=250000/(1+ae^(bt)) and find a and b.         

In 1970, when t=0, y=77435              

Thus, 77435=250000/(1+a)                        
So that 1+a=250,000/77435=                      

Thus, a=2.228514                        

Now, find b.  Use other actual data endpoint from 1987.                 

In 1987, when t=16, y=184790                    

184790=250000/(1+2.228514*exp(16*b))                    
1+2.228514*exp(16*b)=250,000/184,790=                   1.352887061 

exp(16*b)=(1.352887-1)/2.228514=                        0.158350811 

16b=ln(0.158350)=                       0.459637587 

b=(ln(0.158350))/16=                    -0.115184219 


Logistic equation fit using actual data:                        

y=250,000/(1+2.228514e^(-0.11518*t))                    

t       logistic        exponential     actual

0       77435.0057      74332.97938     77435 
1       83722.23968     78484.54458     81483 
2       90252.97004     82867.97851     82077 
3       96995.96015     87496.23125     87569 
4       103915.0238     92382.97616     92759 
5       110969.682      97542.65027     102426 
6       118116.0349     102990.4969     105436 
7       125307.8122     108742.6108     101409 
8       132497.5521     114815.9855     106165 
9       139637.8469     121228.5638     108040 
10      146682.5842     127999.2905     127545 
11      153588.1177     135148.1685     136160 
12      160314.3044     142696.3179     145569 
13      166825.359      150666.0384     148470 
14      173090.491      159080.8751     163950 
15      179084.3105     167965.6881     174324 
16      184787.0014     177346.7261     184790 
17      190184.2805     187251.7037     
18      195267.1697     197709.8834     
19      200031.6182     208752.1621     
20      204478.0126     220411.1622     
21      208610.6158     232721.3281     
22      212436.9698     245719.0281     
23      215967.2948     259442.6616     
24      219213.9089     273932.7727     
25      222190.6871     289232.1697     
26      224912.574      305386.0522     
27      227395.1549     322442.1438     
28      229654.2899     340450.8338     
29      231705.8104     359465.3257     
30      233565.2733     379541.7944     
31      235247.7695     400739.5524     
32      236767.7792     423121.2246     
33      238139.0698     446752.9338     
34      239374.6272     471704.4957     
35      240486.618      498049.6252     
36      241486.3731     525866.1544     
37      242384.3918     555236.2623     
38      243190.3583     586246.7177     
39      243913.1707     618989.1356     
40      244560.9758     653560.2476     
41      245141.2101     690062.188      
42      245660.6434     728602.7952     
43      246125.4239     769295.9308     
44      246541.1236     812261.8155     
45      246912.7826     857627.3842     
46      247244.9529     905526.6617     
47      247541.739      956101.1578     
48      247806.8375     1009500.286     
49      248043.5738     1065881.805     
50      248254.9361     1125412.283     
51      248443.6071     1188267.593     
52      248611.9937     1254633.43      
53      248762.2532     1324705.86      
54      248896.3185     1398691.899     

Suppose instead, that the 1987 data was suspect and that you felt                       
more confident using the 1980 data.                     

Still, a=2.228514                       
But, in 1980--when t=9, y=108040.                       

So, 108,040=250,000/(1+2.228514*e^(9b))                 
(1+2.228512*e^(9b))=250,000/108,040=                    2.313957793 

Thus, e^(9b)=(2.313957-1)/2.228512=                     0.589611813 

So, 9b=ln(0.589611)                     -0.528292282 
b=                      -0.058699142 

So, the logistic equation in this case is:                              
     y=250,000/(1+2.228514e^(-0.05869914*t))                            

t       logistic 2      logistic 1      exponential     actual
0       77435.0057      77435.0057      74332.97938     77435 
1       80607.00596     83722.23968     78484.54458     81483 
2       83845.8089      90252.97004     82867.97851     82077 
3       87147.79645     96995.96015     87496.23125     87569 
4       90508.98645     103915.0238     92382.97616     92759 
5       93925.04538     110969.682      97542.65027     102426 
6       97391.3053      118116.0349     102990.4969     105436 
7       100902.7849     125307.8122     108742.6108     101409 
8       104454.2146     132497.5521     114815.9855     106165 
9       108040.0652     139637.8469     121228.5638     108040 
10      111654.5804     146682.5842     127999.2905     127545 
11      115291.8122     153588.1177     135148.1685     136160 
12      118945.6588     160314.3044     142696.3179     145569 
13      122609.9054     166825.359      150666.0384     148470 
14      126278.2657     173090.491      159080.8751     163950 
15      129944.4254     179084.3105     167965.6881     174324 
16      133602.085      184787.0014     177346.7261     184790 
17      137245.0036     190184.2805     187251.7037     
18      140867.0408     195267.1697     197709.8834     
19      144462.1979     200031.6182     208752.1621     
20      148024.6565     204478.0126     220411.1622     
21      151548.8151     208610.6158     232721.3281     
22      155029.3221     212436.9698     245719.0281     
23      158461.1061     215967.2948     259442.6616     
24      161839.4017     219213.9089     273932.7727     
25      165159.7721     222190.6871     289232.1697     
26      168418.1269     224912.574      305386.0522     
27      171610.7368     227395.1549     322442.1438     
28      174734.2426     229654.2899     340450.8338     
29      177785.662      231705.8104     359465.3257     
30      180762.3905     233565.2733     379541.7944     
31      183662.2007     235247.7695     400739.5524     
32      186483.2364     236767.7792     423121.2246     
33      189224.0046     238139.0698     446752.9338     
34      191883.3646     239374.6272     471704.4957     
35      194460.5141     240486.618      498049.6252     
36      196954.974      241486.3731     525866.1544     
37      199366.571      242384.3918     555236.2623     
38      201695.4187     243190.3583     586246.7177     
39      203941.8984     243913.1707     618989.1356     
40      206106.6384     244560.9758     653560.2476     
41      208190.4932     245141.2101     690062.188      
42      210194.5229     245660.6434     728602.7952     
43      212119.9724     246125.4239     769295.9308     
44      213968.2506     246541.1236     812261.8155     
45      215740.9111     246912.7826     857627.3842     
46      217439.6328     247244.9529     905526.6617     
47      219066.2015     247541.739      956101.1578     
48      220622.4923     247806.8375     1009500.286     
49      222110.4537     248043.5738     1065881.805     
50      223532.0918     248254.9361     1125412.283     
51      224889.4562     248443.6071     1188267.593     
52      226184.6266     248611.9937     1254633.43      
53      227419.7012     248762.2532     1324705.86      
54      228596.7854     248896.3185     1398691.899