使用Python实现SEIR模型,模拟流行病传播过程
1.人群划分:在本模型中,人群四种状态:易感者(S), 暴露者(E),感染者(I)以及康复者(R) 2.模型假设人群分布是同质均匀的,未考虑人群出生、死亡、迁入迁出对疾病传播的影响 3.康复者永久免疫:康复者永久免疫该传染病,无再感染风险 4.干预措施:戴口罩会降低易感者被传染的概率
以20XX年XX月至20XX年XX月在中国各省的某传染病病例负担作为基准,拟合和调整模型参数,构建可用的SEIR模型(易感者-暴露者-感染者康复者/移除者)。并且考虑在干预措施(戴口罩)的影响下,S和E人群状态的变化。
参数 | 定义 | 单位 |
---|---|---|
β | 人传人传播系数 | |
ω | 潜伏期速率 | day-1 |
γ | 恢复期速率 | day-1 |
$$
\begin{aligned}
\frac{dS}{dt} &= -\frac{\beta IS}{N} \\
\frac{dE}{dt} &= \frac{\beta IS}{N} - \omega E \\
\frac{dI}{dt} &= \omega E - \gamma I \\
\frac{dR}{dt} &= \gamma I \\
N &= S + E + I + R
\end{aligned}
$$
适用于疾病病程中存在易感者-暴露者-感染者-康复者/移除者四种状态的传染病。
本SEIR模型模拟了某传染病的流行趋势,并通过添加干预措施(如戴口罩),模拟了戴口罩人群比例为0.00%, 20%, 40%, 60%, 80%的情况下,该传染病的流行趋势,从而更加直观的得到不同程度干预措施下的疾病负担变化情况和干预效果,为政策指定提供参考。
适用于存在疾病病程中存在易感者-暴露者-感染者-康复者/移除者四种状态的传染病。
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp
from lmfit import Model, Parameters
import pandas as pd
# 设置Matplotlib支持中文显示
plt.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['axes.unicode_minus'] = False
# 定义 SEIR 模型的微分方程
def seir_model(t, y, N, beta, omega, gamma):
S, E, I, R = y
dS_dt = -beta * S * I / N
dE_dt = beta * S * I / N - omega * E
dI_dt = omega * E - gamma * I
dR_dt = gamma * I
return [dS_dt, dE_dt, dI_dt, dR_dt]
# 导入病例数据
observed_I = pd.read_excel(r'C:\Users\xieyulun\Desktop\案例汇报\seir\\\\病例数.xlsx')
observed_I = observed_I.reset_index(drop=True)
# 提取 cases 列的数据
observed_I = observed_I['cases'].values
days = len(observed_I) # 确定时间范围
t = np.linspace(0, days - 1, days) # 时间点
# 初始条件和总人口
N = 10000 # 总人口数量
I0 = observed_I[0] # 初始感染者人数
E0 = I0 * 2 # 假设初始潜伏者人数为感染人数的两倍
R0 = 0 # 初始康复者人数
S0 = N - I0 - E0 - R0 # 初始易感者人数
# 初始状态
y0 = [S0, E0, I0, R0]
# 潜伏率和康复率
omega = 1/5 # 潜伏率
gamma = 1/12 # 康复率
# 创建参数对象
params = Parameters()
params.add('beta', min=0.01, max=1.00)
# 创建拟合函数。
def fit_func(t, beta, N, omega, gamma, y0):
# 求解 ODE
result = solve_ivp(seir_model, [t[0], t[-1]], y0, args=(N, beta, omega, gamma), t_eval=t)
return result.y[2]
# 创建拟合模型对象,使用 fit_func 来拟合数据
fit_model = Model(fit_func, independent_vars=['t', 'N', 'omega', 'gamma', 'y0'])
# 拟合模型
result = fit_model.fit(observed_I, params, t=t, N=N, omega=omega, gamma=gamma, y0=y0)
# 打印拟合结果的报告,包括最优参数值、拟合的统计信息等。
print(result.fit_report())
# 提取最优 beta
optimal_beta = result.params['beta'].value
print(f"最优 beta: {optimal_beta:.4f}")
# 模型预测结果
fitted_solution = solve_ivp(seir_model, [t[0], t[-1]], y0, args=(N, optimal_beta, omega, gamma), t_eval=t)
S, E, I, R = fitted_solution.y
# 绘制结果
plt.figure(figsize=(12, 6))
plt.plot(t, S, label='易感者 S', color='blue')
plt.plot(t, E, label='潜伏者 E', color='orange')
plt.plot(t, I, label='感染者 I', color='red')
plt.plot(t, R, label='康复者 R', color='green')
plt.scatter(t, observed_I, label='观测感染数据', edgecolors='black', facecolors='none', linewidth=1.5)
plt.title("SEIR 模型拟合与感染数据")
plt.xlabel("时间 (天)")
plt.ylabel("人数")
plt.legend()
# 获取当前轴对象
ax = plt.gca()
# 去掉上边框和右边框
ax.spines['top'].set_visible(False)
ax.spines['right'].set_visible(False)
plt.show()
[[Model]]
Model(fit_func)
[[Fit Statistics]]
# fitting method = leastsq
# function evals = 54
# data points = 100
# variables = 1
chi-square = 54231245.3
reduced chi-square = 547790.357
Akaike info crit = 1322.35976
Bayesian info crit = 1324.96493
R-squared = 0.77449417
[[Variables]]
beta: 0.44017782 +/- 0.01482789 (3.37%) (init = 0.01)
最优 beta: 0.4402
# 定义 SEIR 模型的微分方程
def seir_model(t, y, N, beta, omega, gamma, mask_usage, mask_effectiveness):
S, E, I, R = y
beta_effective = beta * (1 - mask_usage * mask_effectiveness)
dS_dt = -beta_effective * S * I / N
dE_dt = beta_effective * S * I / N - omega * E
dI_dt = omega * E - gamma * I
dR_dt = gamma * I
return [dS_dt, dE_dt, dI_dt, dR_dt]
# 定义戴口罩的比例
mask_usages = [0.0, 0.2, 0.4, 0.6, 0.8]
mask_effectiveness = 0.5 # 假设戴口罩可以降低 50% 的感染率
# 定义颜色渐变范围
colors = {
'E': plt.cm.Blues,
'I': plt.cm.Reds,
}
# 绘制结果
plt.figure(figsize=(12, 6))
# 为每个戴口罩比例进行预测并绘制结果
for mask_usage in mask_usages:
# 模型预测结果
fitted_solution = solve_ivp(seir_model, [t[0], t[-1]], y0, args=(N, optimal_beta, omega, gamma, mask_usage, mask_effectiveness), t_eval=t)
S, E, I, R = fitted_solution.y
# 选择渐进色
color_E = colors['E'](0.7 + 0.3 * (mask_usage / max(mask_usages)))
color_I = colors['I'](0.7 + 0.3 * (mask_usage / max(mask_usages)))
# 绘制各个状态的结果
plt.plot(t, E, label=f'戴口罩比例 {mask_usage*100}%,潜伏者 E', linestyle='-', color=color_E)
plt.plot(t, I, label=f'戴口罩比例 {mask_usage*100}%,感染者 I', linestyle='-', color=color_I)
plt.title("SEIR 模型预测与感染数据(考虑戴口罩干预)")
plt.xlabel("时间 (天)")
plt.ylabel("人 数")
plt.ylim(0, 5000)
plt.legend()
ax = plt.gca()
ax.spines['top'].set_visible(False)
ax.spines['right'].set_visible(False)
plt.show()
在干预措施(戴口罩)不断加强的影响下,E和I人群的感染人数峰值会降低,也会延缓感染人数峰值的到来。