文章/答案/技术大牛

发布

社区首页 >专栏 >python实现logistic增长模型、多项式模型

python实现logistic增长模型、多项式模型

悟乙己

发布于 2022-11-30 06:45:21

2.3K01

代码可运行

文章被收录于专栏：素质云笔记素质云笔记

运行总次数：1

代码可运行

文章目录

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()

3 案例：疫情数据拟合

案例来源，很好地把一元二次式拟合和一元三次式拟合，还有高斯函数进行拟合： Covid-19-data-fitting-and-prediction

3.1 案例简述

新冠疫情期间，运用 python，基于疫情相关数据设计了几款疫情预测模型，结果曲线能够很好地与国内疫情发展情况拟合并能较好地预测病例增长的拐点时间。由于湖北疑似数据较多，确诊数据准确性较差，我选择了全国除湖北外确诊人数的数据进行拟合，数据来自@人民日报微博每日发布，把1月21日作为统计第一天，进行数据收集。首先，根据国除湖北外确诊人数数据画出了散点图和折线图。

然后，分别利用python的np.polyfit 和 np.polyld分别进行一元二次式拟合和一元三次式拟合，发现一元三次式拟合程度更高。

通过求解系数矩阵，分别计算出一元二次式系数和一元三次式系数。在钟南山院士提出拐点后，尝试预测拐点。选择了高斯函数模型，利用python的curve_fit对每日增长的确诊数量进行拟合，预测拐点。

3.2 高斯函数详细解读

此时案例中的高斯函数代码为：

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
import math


#自定义高斯函数
def func(x, a,u, sig):
    return  a*(np.exp(-(x - u) ** 2 /(2* sig **2))/(math.sqrt(2*math.pi)*sig))


#定义x、y散点坐标
x = [5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29]
x=np.array(x)
print('x is :\n',x)
num = [365,398,480,619,705,761,752,668,722,888,730,707,696,544,505,442,370,377,311,267,221,165,115,81,56]
y = np.array(num)
print('y is :\n',y)

popt, pcov = curve_fit(func, x, y,p0=[0,1,2])
#获取popt拟合系数
a = popt[0]
u = popt[1]
sig = popt[2]

yvals = func(x,a,u,sig) #拟合y值
print('系数a:', a)
print('系数u:', u)
print('系数sig:', sig)
#绘图
plt.rcParams['font.sans-serif']=['SimHei']
plt.xlabel('统计天数')
plt.ylabel('人数')
plot1 = plt.plot(x, y, 's',label='每日增加确诊患者人数',color='blue')
x = np.linspace(u - 3*sig, u + 3*sig, 50)
x_01 = np.linspace(u - 6 * sig, u + 6 * sig, 50)
y_sig = a* ( np.exp(-(x - u) ** 2 /(2* sig **2))/(math.sqrt(2*math.pi)*sig) )
plt.plot(x, y_sig, "r-", linewidth=2)
# plt.plot(x, y, 'r-', x, y, 'go', linewidth=2,markersize=8)
plt.grid(True)
plt.legend(loc=2)
plt.title('每日增长高斯曲线')
plt.show()

图示为：

首先这里的高斯函数：

def func(x, a,u, sig,penalty = 0.8):
    '''
    penalty 惩罚项,越小,尖峰越高
    '''
    return  a*(np.exp(-(x - u) ** 2 /(2* sig **2))/(math.sqrt(2*math.pi)*sig))*penalty

笔者自己改了之后，新增了一个惩罚项，其实就是原来的a，一般来说：

penalty越小，顶峰越尖
penalty越大，顶峰越矮

这个就是惩罚项很小0.1时候的结果

另外，这里生成预测集合的时候：

x = np.linspace(u - 3*sig, u + 3*sig, 50)
x_01 = np.linspace(u - 6 * sig, u + 6 * sig, 50)
y_sig = a* ( np.exp(-(x - u) ** 2 /(2* sig **2))/(math.sqrt(2*math.pi)*sig) )

可以不用这么麻烦，还得用linspace 笔者最终微调之后的代码为：

x = np.array(df_fb['ds'])  #.reshape(-1, 1)# df_fb['ds']
y = np.array(df_fb['y'])#.reshape(-1, 1)

#自定义高斯函数
def func(x, a,u, sig,penalty = 0.8):
    '''
    penalty 惩罚项,越小,尖峰越高
    '''
    return  a*(np.exp(-(x - u) ** 2 /(2* sig **2))/(math.sqrt(2*math.pi)*sig))*penalty

popt, pcov = curve_fit(func, x, y,p0=[0,1,2])

#获取popt拟合系数
a = popt[0]
u = popt[1]
sig = popt[2]


yvals = func(x,a,u,sig) #拟合y值

#绘图
plt.rcParams['font.sans-serif']=['SimHei']
plt.xlabel('统计天数')
plt.ylabel('人数')
plot1 = plt.plot(x, y, 's',label='每日增加确诊患者人数',color='blue')

# 生成预测集合
x_random = list(range(len(data)))
y_sig = a* ( np.exp(-(x_random - u) ** 2 /(2* sig **2))/(math.sqrt(2*math.pi)*sig) )
plt.plot(x_random, y_sig, "r-", linewidth=2)

plt.grid(True)
plt.legend(loc=2)
plt.title('每日增长高斯曲线')
plt.show()