原先写了篇论文,现在把代码分享给大家。
(谢谢大家的关注,由于某些原因,请原谅相关数据信息不能共享;工作繁忙,研究久远,请原谅相关问题无法解答)
代码
import scipy.integrate as spi
import numpy as np
import matplotlib.pyplot as plt
import math
import xlrd
# https://zhuanlan.zhihu.com/p/104091330
# I_0为感染者的初始人数 524*1.5
I_0 = 109
# E_0为潜伏者的初始人数
E_0 = 223
# R_0为治愈者的初始人数
R_0 = 4
# S_0为易感者的初始人数
# S_0 = N - I_0 - E_0 - R_0
S_0=21536000
# Sq 隔离易感染者
S_q=1000
# Eq隔离潜伏者
E_q=200
# H住院患者
H=I_0+E_q
# T为传播时间
T = 35
#θ潜伏者相对感染者传播能力的比值-假设潜伏期患者与已经表现出症状的患者传染能力相同
cita=1
#γ隔离接触速率-隔离14天
gamma=1/14
# σ潜伏者向感染者转化速率 潜伏期7天
sigama=1/7
# α病死率----需基于实际死亡人数进行调整
alpha=0.013
# δI感染者隔离速率
delta_I=0.13
# γI感染者恢复速率
gamma_I=0.043
# δq 隔离潜伏者向隔离感染者的转化速率
delta_q=0.013
# γH隔离感染者的恢复速率
gamma_H=0.043
# β 传染概率
beta = 3.3*math.pow(10,-9)
# q 隔离比例
q=59*math.pow(10,-6)
# c 接触率
c=2.9
# ρ有效接触系数 ρc有效接触率
rou=1
#λ隔离解除速率
lamada=1/14
# INI为初始状态下的数组
INI = (S_0,E_0,I_0,R_0,S_q,E_q,H)
def funcSEIR(inivalue,_):
Y = np.zeros(7)
X = inivalue
Y[0]=-(rou*c*beta+rou*c*q*(1-beta))*X[0]*(X[2]+cita*X[1])+lamada*X[4]
Y[1]=rou*c*beta*(1-q)*X[0]*(X[2]+cita*X[1])-sigama*X[1]
Y[2]=sigama*X[1]-(delta_I+alpha+gamma_I)*X[2]
Y[3]=rou*c*q*(1-beta)*X[0]*(X[2]+cita*X[1])-lamada*X[4]
Y[4]=rou*c*beta*q*X[0]*(X[2]+cita*X[1])-delta_I*X[4]
Y[5]=delta_I*X[2]+delta_q*X[5]-(alpha+gamma_H)*X[6]
Y[6]=gamma_I*X[2]+gamma_H*X[6]
return Y
T_range = np.arange(0,T + 1)
RES = spi.odeint(funcSEIR,INI,T_range)
# plt.plot(RES[:,0],color = 'darkblue',label = 'Susceptible',marker = '.')
# plt.plot(RES[:,1],color = 'orange',label = 'Exposed',marker = '.')
plt.plot(RES[:,2],color = 'red',label = '模拟感染人数',marker = '.')
# plt.plot(RES[:,3],color = 'green',label = 'S_q',marker = '.')
# plt.plot(RES[:,4],color = 'blue',label = 'E_q',marker = '.')
# plt.plot(RES[:,5],color = 'black',label = 'H',marker = '.')
# plt.plot(RES[:,6],color = 'darkgreen',label = 'Recovery',marker = '.')
filename=r'paint1.xlsx'
file= xlrd.open_workbook(filename) # 打开文件
table=file.sheets()[3]#通过索引顺序获取 从0开始
x=table.row_values(2)
plt.rcParams['font.sans-serif']=['SimHei']#用来正常显示中文标签
plt.plot(x,label='现实感染人数')
plt.title('北京')
plt.legend()
plt.xlabel('天数')
plt.ylabel('人数')
plt.show()
效果图
参考文献
1. 代码参考-用Python实现经典传染病模型
2. 理论参考-曹盛力,冯沛华,时朋朋.修正SEIR传染病动力学模型应用于湖北省2019冠状病毒病(COVID-19)疫情预测和评估[J/OL].浙江大学学报(医学版):1-13[2020-03-29].http://kns.cnki.net/kcms/detail/33.1248.R.20200303.1722.004.html.