用Python从零复现传染病模型SI、SIS、SIR、SEIRD的保姆级代码实现与对比传染病模型是理解疾病传播动态的重要工具。对于开发者、数据分析师和学生来说将这些数学模型转化为可运行的代码不仅能加深理解还能为实际应用提供有力支持。本文将带你用Python的NumPy、SciPy和Matplotlib库从零开始实现四种经典传染病模型SI、SIS、SIR和SEIRD并通过可视化对比不同参数下的传播曲线。1. 环境准备与基础概念在开始编码前我们需要确保环境配置正确并理解基本概念。Python的科学计算生态系统为我们提供了强大的工具集# 必需库安装如未安装 # pip install numpy scipy matplotlib pandas传染病模型的核心是微分方程描述人群在不同状态间的转移。基本参数包括传染率(λ)一个感染者每天能传染多少易感者治愈率(μ)感染者每天被治愈的概率基本再生数(R₀)一个感染者在易感人群中能传染的平均人数提示在Jupyter Notebook中运行代码时建议先导入所有需要的库保持环境整洁。import numpy as np from scipy.integrate import odeint import matplotlib.pyplot as plt plt.style.use(ggplot) # 使用更美观的绘图样式2. SI模型实现最简单的传染病传播SI模型将人群分为两类易感者(S)和感染者(I)。假设疾病不可治愈所有易感者最终都会被感染。2.1 数学模型建立SI模型的微分方程为dS/dt -λSI/N dI/dt λSI/N其中N是总人口数(SIN)。def si_model(y, t, N, lambda_): S, I y dSdt -lambda_ * S * I / N dIdt lambda_ * S * I / N return [dSdt, dIdt]2.2 参数设置与模拟让我们设置初始条件和参数进行模拟# 参数设置 N 1000 # 总人口 I0 1 # 初始感染者 S0 N - I0 # 初始易感者 lambda_ 0.3 # 传染率 t np.linspace(0, 50, 50) # 50天的时间范围 # 解微分方程 solution odeint(si_model, [S0, I0], t, args(N, lambda_)) S, I solution.T2.3 结果可视化plt.figure(figsize(10,6)) plt.plot(t, S, b, label易感者) plt.plot(t, I, r, label感染者) plt.xlabel(时间(天)) plt.ylabel(人数) plt.title(SI模型传播动态) plt.legend() plt.grid(True) plt.show()注意SI模型适合描述不可治愈的传染病如某些慢性病毒感染。在实际应用中传染率λ的估计对模型准确性至关重要。3. SIS模型考虑治愈但无免疫力的情形SIS模型在SI基础上增加了治愈过程但治愈者不获得免疫力可能再次被感染。3.1 模型建立与实现SIS模型的微分方程为dS/dt -λSI/N μI dI/dt λSI/N - μIdef sis_model(y, t, N, lambda_, mu): S, I y dSdt -lambda_ * S * I / N mu * I dIdt lambda_ * S * I / N - mu * I return [dSdt, dIdt]3.2 关键参数对比我们比较不同治愈率下的传播动态参数组合λ值μ值R₀(λ/μ)预期结果情形10.40.14疫情持续情形20.20.30.67疫情消退# 模拟两种情形 mu1, mu2 0.1, 0.3 sol1 odeint(sis_model, [S0, I0], t, args(N, lambda_, mu1)) sol2 odeint(sis_model, [S0, I0], t, args(N, lambda_, mu2))3.3 可视化对比plt.figure(figsize(12,8)) plt.subplot(2,1,1) plt.plot(t, sol1[:,1], r, labelf感染者(R₀{lambda_/mu1:.2f})) plt.title(SIS模型R₀1情形) plt.legend() plt.subplot(2,1,2) plt.plot(t, sol2[:,1], g, labelf感染者(R₀{lambda_/mu2:.2f})) plt.title(SIS模型R₀1情形) plt.legend() plt.tight_layout() plt.show()4. SIR模型经典传染病框架SIR模型引入了康复者(R)群体他们获得免疫力不再被感染。4.1 模型实现SIR模型的微分方程为dS/dt -λSI/N dI/dt λSI/N - μI dR/dt μIdef sir_model(y, t, N, lambda_, mu): S, I, R y dSdt -lambda_ * S * I / N dIdt lambda_ * S * I / N - mu * I dRdt mu * I return [dSdt, dIdt, dRdt]4.2 参数敏感性分析改变R₀值观察疫情发展R0_values [1.5, 2.5, 3.5] mu 0.2 plt.figure(figsize(12,8)) for R0 in R0_values: lambda_ R0 * mu solution odeint(sir_model, [S0, I0, 0], t, args(N, lambda_, mu)) plt.plot(t, solution[:,1], labelfR₀{R0}) plt.title(不同R₀值下的感染曲线) plt.xlabel(时间(天)) plt.ylabel(感染人数) plt.legend() plt.grid(True) plt.show()4.3 群体免疫阈值SIR模型的一个重要结论是群体免疫阈值HIT 1 - 1/R₀def herd_immunity_threshold(R0): return 1 - 1/R0 R0_range np.linspace(1.1, 10, 100) hit_values [herd_immunity_threshold(r) for r in R0_range] plt.figure(figsize(10,6)) plt.plot(R0_range, hit_values) plt.xlabel(基本再生数(R₀)) plt.ylabel(群体免疫阈值) plt.title(R₀与群体免疫阈值关系) plt.grid(True) plt.show()5. SEIRD模型新冠疫情的复杂模拟SEIRD模型在SIR基础上增加了潜伏期(E)和死亡(D)两个状态更适合COVID-19等现代传染病。5.1 模型建立SEIRD模型的微分方程为dS/dt -βSI/N dE/dt βSI/N - σE dI/dt σE - (γ α)I dR/dt γI dD/dt αIdef seird_model(y, t, N, beta, sigma, gamma, alpha): S, E, I, R, D y dSdt -beta * S * I / N dEdt beta * S * I / N - sigma * E dIdt sigma * E - (gamma alpha) * I dRdt gamma * I dDdt alpha * I return [dSdt, dEdt, dIdt, dRdt, dDdt]5.2 参数估计与设置根据COVID-19研究设置典型参数# 参数设置 N 10000 # 总人口 E0, I0 10, 5 # 初始潜伏者和感染者 R0, D0 0, 0 # 初始康复者和死亡者 S0 N - E0 - I0 - R0 - D0 # 流行病学参数 R0_value 2.5 # 基本再生数 sigma 1/5.2 # 潜伏期倒数(平均5.2天) gamma 1/14 # 治愈率(平均14天) alpha 0.01 # 死亡率(1%) beta R0_value * gamma # 传染率 t np.linspace(0, 180, 180) # 模拟180天5.3 完整模拟与可视化solution odeint(seird_model, [S0, E0, I0, R0, D0], t, args(N, beta, sigma, gamma, alpha)) S, E, I, R, D solution.T plt.figure(figsize(12,8)) plt.plot(t, S, b, label易感者) plt.plot(t, E, y, label潜伏者) plt.plot(t, I, r, label感染者) plt.plot(t, R, g, label康复者) plt.plot(t, D, k, label死亡者) plt.xlabel(时间(天)) plt.ylabel(人数) plt.title(SEIRD模型模拟COVID-19传播) plt.legend() plt.grid(True) plt.show()5.4 关键指标分析计算几个重要指标peak_day t[np.argmax(I)] peak_infections np.max(I) total_deaths D[-1] attack_rate (R[-1] D[-1]) / N * 100 print(f疫情高峰出现在第{peak_day:.0f}天) print(f最高同时感染人数{peak_infections:.0f}) print(f最终死亡人数{total_deaths:.0f}) print(f总感染率{attack_rate:.1f}%)6. 模型对比与应用建议四种模型各有适用场景选择取决于具体需求模型适用场景优点局限性SI不可治愈疾病简单直观过于简化SIS可重复感染疾病考虑治愈忽略免疫力SIR获得免疫力的疾病经典框架忽略潜伏期SEIRD现代复杂传染病全面考虑多个状态参数较多在实际项目中我通常会先使用SIR模型进行快速评估当需要更精确的结果时再转向SEIRD模型。参数估计是模型准确性的关键建议从文献中获取初步参数范围使用实际数据进行参数校准进行敏感性分析确定关键参数用交叉验证评估模型预测能力# 模型对比可视化函数 def compare_models(t, N, lambda_0.3, mu0.1, sigma1/5, gamma1/14, alpha0.01): # SI模型 si_sol odeint(si_model, [N-1, 1], t, args(N, lambda_)) # SIS模型 sis_sol odeint(sis_model, [N-1, 1], t, args(N, lambda_, mu)) # SIR模型 sir_sol odeint(sir_model, [N-1, 1, 0], t, args(N, lambda_, mu)) # SEIRD模型 beta lambda_ seird_sol odeint(seird_model, [N-1, 0, 1, 0, 0], t, args(N, beta, sigma, gamma, alpha)) plt.figure(figsize(12,10)) plt.plot(t, si_sol[:,1], r--, labelSI模型感染者) plt.plot(t, sis_sol[:,1], g-., labelSIS模型感染者) plt.plot(t, sir_sol[:,1], b:, labelSIR模型感染者) plt.plot(t, seird_sol[:,2], k-, labelSEIRD模型感染者) plt.title(四种传染病模型感染曲线对比) plt.xlabel(时间(天)) plt.ylabel(感染人数) plt.legend() plt.grid(True) plt.show()在公共卫生决策中这些模型可以帮助评估不同干预措施的效果。例如通过调整传染率参数来模拟社交距离措施的影响或改变治愈率来模拟医疗资源投入的效果。