seir模型数学建模python_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()

效果图

参考文献

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.

  • 0
    点赞
  • 12
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
下面是使用 Gillespie Algorithm 实现 SEIR 模型Python 代码: ```python import numpy as np import matplotlib.pyplot as plt def gillespieSEIR(beta, sigma, gamma, N, E0, I0, R0, Tmax): # Initialize variables t = 0 S, E, I, R = N - E0 - I0 - R0, E0, I0, R0 S_list, E_list, I_list, R_list, t_list = [S], [E], [I], [R], [t] while t < Tmax: # Calculate rates of each event lambda1 = beta*S*I/N lambda2 = sigma*E lambda3 = gamma*I # Calculate total rate lambda_total = lambda1 + lambda2 + lambda3 # Sample time until next event dt = np.random.exponential(scale=1/lambda_total) # Sample which event occurs event = np.random.choice(['lambda1', 'lambda2', 'lambda3'], p=[lambda1/lambda_total, lambda2/lambda_total, lambda3/lambda_total]) # Update variables if event == 'lambda1': S -= 1 I += 1 elif event == 'lambda2': E -= 1 I += 1 else: I -= 1 R += 1 # Update time and lists t += dt S_list.append(S) E_list.append(E) I_list.append(I) R_list.append(R) t_list.append(t) return S_list, E_list, I_list, R_list, t_list if __name__ == '__main__': # Set parameters beta = 0.2 sigma = 0.1 gamma = 0.05 N = 1000 E0, I0, R0 = 10, 1, 0 Tmax = 100 # Run simulation S, E, I, R, t = gillespieSEIR(beta, sigma, gamma, N, E0, I0, R0, Tmax) # Plot results plt.plot(t, S, label='Susceptible') plt.plot(t, E, label='Exposed') plt.plot(t, I, label='Infected') plt.plot(t, R, label='Recovered') plt.legend() plt.xlabel('Time') plt.ylabel('Number of individuals') plt.show() ``` 其中,参数 `beta` 表示感染率,`sigma` 表示潜伏期转化为感染期的速率,`gamma` 表示康复率,`N` 表示总人口数,`E0`、`I0`、`R0` 表示初始时刻的潜伏期、感染期和康复人数,`Tmax` 表示模拟的时间长度。函数 `gillespieSEIR` 返回每个时间点的 S、E、I、R 数量,以及对应的时间点。最后,使用 matplotlib 库将结果可视化。

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值