1 logistic growth model
1.1 J-type growth and S-type growth
Exponential growth, J-shaped curve: exponential growth, that is, the growth is not restrained and explosive.
For example, one person can infect three people, three people can infect nine people, and nine people can infect 27 people, which keeps doubling. This is J-type growth, also known as exponential growth.
Some infectious diseases may show exponential growth in the initial stage.
However, in the actual growth process, the growth rate cannot be maintained all the time. With the increasing number of people, the growth rate will be gradually restrained. This is S-type growth.
The spread of general diseases is an S-type growth process, because there will be some resistance in the process of disease transmission.
1.2 logistic growth function
When a species moves into a new ecosystem, its number will change. Assuming that the initial number of the species is less than the maximum capacity of the environment, the number will increase. If the species has natural enemies, food, space and other resources in this ecosystem and is also insufficient (non ideal environment), the growth function meets the logistic equation, and the image is S-shaped. This equation is the best mathematical model to describe the law of population growth under the condition of limited resources. In the following contents, the principle, ecological significance and application of logistic equation will be introduced in detail. The differential formula of logistic model is: dx/dt=rx(1-x), in which r is the rate parameter.
- K is the environmental capacity, that is, the limit that P(t) can reach when it grows to the end.
- P0 is the initial capacity, which is the number at t=0.
- r is the growth rate. The larger r is, the faster it will grow and approach the value of K. The smaller r is, the slower it will grow and approach the value of K.
The formula is written in python in the form of function:
def logistic_increase_function(t,K,P0,r):
# t:time t0:initial time P0:initial_value K:capacity r:increase_rate
exp_value=np.exp(r*(t-t0))
return (K*exp_value*P0)/(K+(exp_value-1)*P0)
The curve of logistic growth is also called s-shaped curve. The figure below shows the number of curves on the left and the growth rate on the right.
1.3 case code
#!/usr/bin/python
# -*- coding: UTF-8 -*-
"""
Fitting 2019-nCov Number of confirmed pneumonia infections
"""
import numpy as np
import matplotlib.pyplot as plt
import math
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
def logistic_increase_function(t,K,P0,r):
'''
t ,list,Date series,[11,18,19,20 ,21, 22, 23, 24, 25, 26, 27]
t0,int,Date first day
r,float,r Is the growth rate, r The larger, the faster the growth, the faster the approach K value - Steeper;r The smaller, the slower the growth, the slower the approach K Value.
P0 Is the initial capacity, which is t=0 Number of moments
K,float,K Is the environmental capacity, that is, to grow to the end, P(t)The limit that can be reached,Generally 1
'''
t0=11 # first day
r=0.6
# r = 0.55
# t:time t0:initial time P0:initial_value K:capacity r:increase_rate
exp_value=np.exp(r*(t-t0))
return (K*exp_value*P0)/(K+(exp_value-1)*P0)
'''
1.11 41 cases per day
1.18 45 cases per day
1.19 62 cases per day
1.20 291 cases per day
1.21 440 cases per day
1.22 571 cases per day
1.23 830 cases per day
1.24 1287 cases per day
1.25 1975 cases in Japan
1.26 2744 cases per day
1.27 4515 cases per day
'''
# Date and number of infected persons
t=[11,18,19,20 ,21, 22, 23, 24, 25, 26, 27]
t=np.array(t)
P=[41,45,62,291,440,571,830,1287,1975,2744,4515]
P=np.array(P)
# Estimation and fitting by least square method
popt, pcov = curve_fit(logistic_increase_function, t, P)
# popt - K，P0，r
# Finally, K=4.01665705e+10 people will be infected
# array([7.86278276e+03, 2.96673434e-01, 1.00000000e+00])
#Get the fitting coefficient in popt
print("K:capacity P0:initial_value r:increase_rate t:time")
print(popt)
#Predicted P value after fitting
P_predict = logistic_increase_function(t,popt[0],popt[1],popt[2])
#Future forecast
future=[11,18,19,20 ,21, 22, 23, 24, 25, 26, 27,28,29,30,31,41,51,61,71,81,91,101]
future=np.array(future)
future_predict=logistic_increase_function(future,popt[0],popt[1],popt[2])
#Recent situation forecast
tomorrow=[28,29,30,32,33,35,37,40]
tomorrow=np.array(tomorrow)
tomorrow_predict=logistic_increase_function(tomorrow,popt[0],popt[1],popt[2])
#mapping
plot1 = plt.plot(t, P, 's',label="confimed infected people number")
plot2 = plt.plot(t, P_predict, 'r',label='predict infected people number')
plot3 = plt.plot(tomorrow, tomorrow_predict, 's',label='predict infected people number')
plt.xlabel('time')
plt.ylabel('confimed infected people number')
plt.legend(loc=0) #Specify the lower right corner of the position of legend
print(logistic_increase_function(np.array(28),popt[0],popt[1],popt[2]))
print(logistic_increase_function(np.array(29),popt[0],popt[1],popt[2]))
plt.show()
#Future forecast mapping
#plot2 = plt.plot(t, P_predict, 'r',label='polyfit values')
#plot3 = plt.plot(future, future_predict, 'r',label='polyfit values')
#plt.show()
print("Program done!")
'''
# Fitting age
'''
import numpy as np
import matplotlib.pyplot as plt
#Define x and y scatter coordinates
x = [10,20,30,40,50,60,70,80]
x = np.array(x)
print('x is :\n',x)
num = [174,236,305,334,349,351,342,323]
y = np.array(num)
print('y is :\n',y)
#Fitting with cubic polynomial
f1 = np.polyfit(x, y, 3)
print('f1 is :\n',f1)
p1 = np.poly1d(f1)
print('p1 is :\n',p1)
#You can also use yvals = NP polyval(f1, x)
yvals = p1(x)
print('yvals is :\n',yvals)
#mapping
plot1 = plt.plot(x, y, 's',label='original values')
plot2 = plt.plot(x, yvals, 'r',label='polyfit values')
plt.xlabel('x')
plt.ylabel('y')
plt.legend(loc=4) #Specify the position of the lower right corner of the legend
plt.title('polyfitting')
plt.show()
The least square method is used for fitting, and the least square method is used to fit the logistic growth function.
logistic_ increase_ The value of r in function (T, K, P0, r) can be adjusted:
After human intervention, the disease reduces the K value, so the r value can be increased to speed up the speed of reaching the K value
(r becomes larger and the curve becomes steeper)
r = 0.55
r=0.65
2 fitting polynomial function
reference resources: Three schemes of quasi merging arbitrary data and curves to obtain function expression in python.
2.1 polynomial fitting - polyfit fitting age
import numpy as np
import matplotlib.pyplot as plt
#Define x and y scatter coordinates
x = [10,20,30,40,50,60,70,80]
x = np.array(x)
print('x is :\n',x)
num = [174,236,305,334,349,351,342,323]
y = np.array(num)
print('y is :\n',y)
#Fitting with cubic polynomial
f1 = np.polyfit(x, y, 3)
print('f1 is :\n',f1)
p1 = np.poly1d(f1)
print('p1 is :\n',p1)
#You can also use yvals = NP polyval(f1, x)
yvals = p1(x) #Fitting y value
print('yvals is :\n',yvals)
#mapping
plot1 = plt.plot(x, y, 's',label='original values')
plot2 = plt.plot(x, yvals, 'r',label='polyfit values')
plt.xlabel('x')
plt.ylabel('y')
plt.legend(loc=4) #Specify the lower right corner of the position of legend
plt.title('polyfitting')
plt.show()
2.2 polynomial fitting - curve_fit polynomial
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
#Custom function e exponential form
def func(x, a, b,c):
return a*np.sqrt(x)*(b*np.square(x)+c)
#Define x and y scatter coordinates
x = [20,30,40,50,60,70]
x = np.array(x)
num = [453,482,503,508,498,479]
y = np.array(num)
#Nonlinear least square fitting
popt, pcov = curve_fit(func, x, y)
#Get the fitting coefficient in popt
print(popt)
a = popt[0]
b = popt[1]
c = popt[2]
yvals = func(x,a,b,c) #y value fitting
print('popt:', popt)
print('coefficient a:', a)
print('coefficient b:', b)
print('coefficient c:', c)
print('coefficient pcov:', pcov)
print('coefficient yvals:', yvals)
#mapping
plot1 = plt.plot(x, y, 's',label='original values')
plot2 = plt.plot(x, yvals, 'r',label='polyfit values')
plt.xlabel('x')
plt.ylabel('y')
plt.legend(loc=4) #Specify the lower right corner of the position of legend
plt.title('curve_fit')
plt.show()
2.3 curve_fit fitting Gaussian distribution
#encoding=utf-8
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
import pandas as pd
#Custom function e exponential form
def func(x, a,u, sig):
return a*(np.exp(-(x - u) ** 2 /(2* sig **2))/(math.sqrt(2*math.pi)*sig))*(431+(4750/x))
#Define x and y scatter coordinates
x = [40,45,50,55,60,65,70,75,80,85,90,95,100,105,110,115,120,125,130,135]
x=np.array(x)
# x = np.array(range(20))
print('x is :\n',x)
num = [536,529,522,516,511,506,502,498,494,490,487,484,481,478,475,472,470,467,465,463]
y = np.array(num)
print('y is :\n',y)
popt, pcov = curve_fit(func, x, y,p0=[3.1,4.2,3.3])
#Get the fitting coefficient in popt
a = popt[0]
u = popt[1]
sig = popt[2]
yvals = func(x,a,u,sig) #Fitting y value
print(u'coefficient a:', a)
print(u'coefficient u:', u)
print(u'coefficient sig:', sig)
#mapping
plot1 = plt.plot(x, y, 's',label='original values')
plot2 = plt.plot(x, yvals, 'r',label='polyfit values')
plt.xlabel('x')
plt.ylabel('y')
plt.legend(loc=4) #Specify the lower right corner of the position of legend
plt.title('curve_fit')
plt.show()
