python实现logistic增长模型

Posted 2021-05-19 悟乙己

文章目录

• 1 logistic 增长模型
• • 1.1 J型增长和S型增长
  • 1.2 logistic增长函数
  • 1.3 案例代码
• 2 拟合多项式函数
• • 2.1 多项式拟合 —— polyfit 拟合年龄
  • 2.2 多项式拟合 —— curve_fit拟合多项式
  • 2.3 curve_fit拟合高斯分布

---

1 logistic 增长模型

1.1 J型增长和S型增长

指数增长，J型曲线：指数增长，即增长不受抑制，呈爆炸式的。

比如一个人可以传染三个人，三个人传染九个人，九个人传染27个人，不停的倍增。这就是J型增长，也叫指数型的增长。

一些传染病初期可能呈现指数增长。

但是实际的增长过程中，增长速率并不能一直维持不变，随着人数的不断增多，增长率会逐渐受到抑制。这就是S型增长。

一般疾病的传播是S型增长的过程，因为疾病传播的过程中会受到一定的阻力。

1.2 logistic增长函数

当一个物种迁入到一个新生态系统中后，其数量会发生变化。假设该物种的起始数量小于环境的最大容纳量，则数量会增长。该物种在此生态系统中有天敌、食物、空间等资源也不足（非理想环境），则增长函数满足逻辑斯谛方程，图像呈S形，此方程是描述在资源有限的条件下种群增长规律的一个最佳数学模型。在以下内容中将具体介绍逻辑斯谛方程的原理、生态学意义及其应用。逻辑斯蒂模型的微分式是：dx/dt=rx(1-x) 式中的r为速率参数。

• K为环境容量，即增长到最后，P(t)能达到的极限。
• P0为初始容量，就是t=0时刻的数量。
• r为增长速率，r越大则增长越快，越快逼近K值，r越小增长越慢，越慢逼近K值。
  该公式用python写成函数形式就是：

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)

logistic增长的曲线也称为s型曲线。下图左图为曲线数量，右图为增长速率。

1.3 案例代码

#!/usr/bin/python
# -*- coding: UTF-8 -*-
"""
拟合2019-nCov肺炎感染确诊人数
"""
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,日期序列,[11,18,19,20 ,21, 22, 23, 24,  25,  26,  27]
    t0,int,日期首日
    r,float,r为增长速率，r越大则增长越快，越快逼近K值 - 越陡峭;r越小增长越慢，越慢逼近K值。
    P0为初始容量，就是t=0时刻的数量
    K,float,K为环境容量，即增长到最后，P(t)能达到的极限,一般为1
    
    '''
    t0=11  # 第一天
    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例
1.18日45例
1.19日62例
1.20日291例
1.21日440例
1.22日571例
1.23日830例
1.24日1287例
1.25日1975例
1.26日2744例
1.27日4515例
'''

#  日期及感染人数
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)
 
# 用最小二乘法估计拟合
popt, pcov = curve_fit(logistic_increase_function, t, P)
# popt - K，P0，r
# 最终K=4.01665705e+10人会被感染
# array([7.86278276e+03, 2.96673434e-01, 1.00000000e+00])

#获取popt里面是拟合系数
print("K:capacity  P0:initial_value   r:increase_rate   t:time")
print(popt)


#拟合后预测的P值
P_predict = logistic_increase_function(t,popt[0],popt[1],popt[2])


#未来预测
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])


#近期情况预测
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])
 
#绘图
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) #指定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()
 
 
#未来预测绘图
#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!")


'''
# 拟合年龄
'''
import numpy as np
import matplotlib.pyplot as plt

#定义x、y散点坐标
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)
#用3次多项式拟合
f1 = np.polyfit(x, y, 3)
print('f1 is :\\n',f1)

p1 = np.poly1d(f1)
print('p1 is :\\n',p1)

#也可使用yvals=np.polyval(f1, x)

yvals = p1(x)
print('yvals is :\\n',yvals)


#绘图
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) #指定legend的位置右下角
plt.title('polyfitting')
plt.show()

用最小二乘法进行拟合,最小二乘法，对logistic增长函数进行拟合。

logistic_increase_function(t,K,P0,r)中的r取值是可以调整的：
人为干预后，疾病降低K值，因此可以将r值提升，以加快达到K值的速度
（r变大，曲线变陡峭）

r取0.55

r=0.65

---

2 拟合多项式函数

参考：python 对于任意数据和曲线进行拟合并求出函数表达式的三种方案。

2.1 多项式拟合 —— polyfit 拟合年龄

import numpy as np
import matplotlib.pyplot as plt
 
#定义x、y散点坐标
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)
#用3次多项式拟合
f1 = np.polyfit(x, y, 3)
print('f1 is :\\n',f1)
 
p1 = np.poly1d(f1)
print('p1 is :\\n',p1)
 
#也可使用yvals=np.polyval(f1, x)
yvals = p1(x)  #拟合y值
print('yvals is :\\n',yvals)
#绘图
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) #指定legend的位置右下角
plt.title('polyfitting')
plt.show()

2.2 多项式拟合 —— curve_fit拟合多项式

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
 
#自定义函数 e指数形式
def func(x, a, b,c):
    return a*np.sqrt(x)*(b*np.square(x)+c)
 
#定义x、y散点坐标
x = [20,30,40,50,60,70]
x = np.array(x)
num = [453,482,503,508,498,479]
y = np.array(num)
 
#非线性最小二乘法拟合
popt, pcov = curve_fit(func, x, y)
#获取popt里面是拟合系数
print(popt)
a = popt[0] 
b = popt[1]
c = popt[2]
yvals = func(x,a,b,c) #拟合y值
print('popt:', popt)
print('系数a:', a)
print('系数b:', b)
print('系数c:', c)
print('系数pcov:', pcov)
print('系数yvals:', yvals)
#绘图
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) #指定legend的位置右下角
plt.title('curve_fit')
plt.show()

2.3 curve_fit拟合高斯分布

#encoding=utf-8  
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
import pandas as pd

#自定义函数 e指数形式
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))


#定义x、y散点坐标
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])
#获取popt里面是拟合系数
a = popt[0]
u = popt[1]
sig = popt[2]


yvals = func(x,a,u,sig) #拟合y值
print(u'系数a:', a)
print(u'系数u:', u)
print(u'系数sig:', sig)

#绘图
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) #指定legend的位置右下角
plt.title('curve_fit')
plt.show()