Python-改进SEIR模型

原先写了篇论文，现在把代码分享给大家。

（谢谢大家的关注，由于某些原因，请原谅相关数据信息不能共享;工作繁忙，研究久远，请原谅相关问题无法解答）

代码

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.