别再只会用SIR模型了！从零到一，用Python+Scipy搞定传染病预测（附SEIR模型代码）

从微分方程到疫情预测Python实现SEIR模型的实战指南当COVID-19席卷全球时许多数据科学家发现自己在凌晨三点盯着电脑屏幕试图理解那些看似简单的微分方程如何转化为真实的疫情曲线。这就是传染病模型的魔力——它们用数学语言讲述着病毒传播的故事。本文将带你走进这个数字化的流行病学世界用Python代码将这些抽象方程转化为可视化的预测工具。1. 传染病模型基础与Python环境准备传染病模型本质上是一组描述人群状态变化的微分方程。想象一下城市中的每个人就像棋盘上的棋子而病毒则是改变棋子颜色的规则。SIR模型将人群分为三类易感者(S)、感染者(I)和康复者(R)而SEIR模型则增加了一个暴露者(E)的类别更贴近具有潜伏期的传染病。1.1 必备工具安装开始之前确保你的Python环境已装备以下武器库pip install numpy scipy matplotlib pandas核心工具说明NumPy处理数组运算的瑞士军刀SciPy科学计算的万能工具箱特别是其中的odeint函数Matplotlib数据可视化的画布Pandas处理真实疫情数据的得力助手提示推荐使用Jupyter Notebook进行交互式开发可以实时观察曲线变化1.2 微分方程求解原理传染病模型的核心是常微分方程组(ODEs)。以经典的SIR模型为例dS/dt -βSI/N dI/dt βSI/N - γI dR/dt γI其中β是感染率γ是康复率。Python中我们用scipy.integrate.odeint来求解这类方程组它采用数值方法逐步计算每个时间点的状态值。2. 从SIR到SEIR模型演进与实现2.1 SIR模型基础实现让我们先用Python实现一个标准的SIR模型import numpy as np from scipy.integrate import odeint import matplotlib.pyplot as plt def sir_model(y, t, beta, gamma): S, I, R y dSdt -beta * S * I dIdt beta * S * I - gamma * I dRdt gamma * I return [dSdt, dIdt, dRdt] # 参数设置 N 1000 # 总人口 beta 0.3 # 感染率 gamma 0.1 # 康复率 I0 1 # 初始感染者 S0 N - I0 # 初始易感者 R0 0 # 初始康复者 # 时间点(天) t np.linspace(0, 160, 160) # 求解微分方程 solution odeint(sir_model, [S0, I0, R0], t, args(beta, gamma)) S, I, R solution.T # 可视化 plt.figure(figsize(10,6)) plt.plot(t, S, label易感者) plt.plot(t, I, label感染者) plt.plot(t, R, label康复者) plt.xlabel(天数) plt.ylabel(人数) plt.legend() plt.show()运行这段代码你会看到经典的SIR曲线感染者数量先上升后下降最终大部分人群进入康复状态。2.2 SEIR模型进阶实现对于像COVID-19这样有潜伏期的疾病SEIR模型更为合适。它在SIR基础上增加了暴露者(E)类别def seir_model(y, t, beta, sigma, gamma): S, E, I, R y dSdt -beta * S * I dEdt beta * S * I - sigma * E dIdt sigma * E - gamma * I dRdt gamma * I return [dSdt, dEdt, dIdt, dRdt] # 新增参数 sigma 1/5.2 # 潜伏期倒数(COVID-19平均潜伏期5.2天) E0 0 # 初始暴露者 solution odeint(seir_model, [S0, E0, I0, R0], t, args(beta, sigma, gamma)) S, E, I, R solution.T # 可视化时增加E曲线参数选择技巧β值通常在0.1-0.5之间取决于疾病传染性γ值约为1/病程天数(如COVID-19约为1/14)σ值约为1/平均潜伏期天数3. 模型调参与实际数据拟合3.1 参数敏感性分析传染病模型对参数非常敏感。下表展示了不同β值对峰值感染人数的影响β值峰值感染人数达到峰值天数最终感染比例0.21504580%0.33003090%0.44002295%进行敏感性分析的代码片段for beta in [0.2, 0.3, 0.4]: solution odeint(seir_model, [S0, E0, I0, R0], t, args(beta, sigma, gamma)) I solution[:, 2] peak np.max(I) print(fβ{beta}: 峰值感染人数{peak:.0f})3.2 使用真实数据拟合参数假设我们有某地区的每日新增病例数据real_cases可以通过最小二乘法拟合最优参数from scipy.optimize import minimize def error_function(params): beta, gamma params solution odeint(sir_model, [S0, I0, R0], t, args(beta, gamma)) predicted np.diff(solution[:,1], prependI0) # 计算每日新增 return np.sum((predicted - real_cases)**2) initial_guess [0.3, 0.1] bounds [(0.01, 0.5), (0.01, 0.3)] result minimize(error_function, initial_guess, boundsbounds) optimal_beta, optimal_gamma result.x注意实际应用中需要考虑数据报告延迟、检测能力等因素可能需要更复杂的误差函数4. 模型扩展与应用场景4.1 考虑人口动态变化基本SEIR模型假设人口恒定我们可以引入出生率ν和死亡率μdef seir_dynamic_model(y, t, beta, sigma, gamma, nu, mu): S, E, I, R y N S E I R dSdt nu * N - mu * S - beta * S * I / N dEdt beta * S * I / N - (mu sigma) * E dIdt sigma * E - (mu gamma) * I dRdt gamma * I - mu * R return [dSdt, dEdt, dIdt, dRdt]4.2 空间异质性建模对于大城市疫情模拟可以考虑将城市划分为多个区域每个区域运行SEIR模型并添加区域间迁移项def multi_seir_model(y, t, beta, sigma, gamma, mobility): # y包含所有区域的S,E,I,R状态 n_regions len(y) // 4 derivatives [] for i in range(n_regions): S, E, I, R y[i*4], y[i*41], y[i*42], y[i*43] N S E I R # 计算迁入迁出 migration_S sum(mobility[j][i] * y[j*4] - mobility[i][j] * S for j in range(n_regions) if j ! i) # 类似计算其他状态的迁移... dSdt -beta[i] * S * I / N migration_S # 其他方程... derivatives.extend([dSdt, dEdt, dIdt, dRdt]) return derivatives4.3 干预措施建模社交隔离措施可以表示为随时间变化的β(t)def beta_with_intervention(t): baseline 0.3 if 30 t 60: # 隔离期 return baseline * 0.3 # 减少70%接触 elif t 60: # 逐步放开 return baseline * min(0.3 (t-60)/100, 1) else: return baseline # 在模型中改用beta_with_intervention(t)代替固定beta5. 可视化与结果解读技巧5.1 动态可视化使用Matplotlib的动画功能展示疫情发展from matplotlib.animation import FuncAnimation fig, ax plt.subplots(figsize(10,6)) line_S, ax.plot([], [], label易感者) line_I, ax.plot([], [], label感染者) def init(): ax.set_xlim(0, len(t)) ax.set_ylim(0, N) return line_S, line_I def update(frame): line_S.set_data(t[:frame], S[:frame]) line_I.set_data(t[:frame], I[:frame]) return line_S, line_I ani FuncAnimation(fig, update, frameslen(t), init_funcinit, blitTrue) plt.legend() plt.show()5.2 关键指标计算从模型结果中提取有用指标peak_day t[np.argmax(I)] # 峰值日期 peak_cases np.max(I) # 峰值感染人数 total_infected N - S[-1] # 总感染人数 # 计算基本再生数R0 R0 beta / gamma print(f基本再生数R0 {R0:.2f}) # 计算群体免疫阈值 herd_immunity_threshold 1 - 1/R0 print(f群体免疫阈值 {herd_immunity_threshold:.1%})结果解读要点R0 1表示疫情会扩散达到群体免疫阈值后疫情自然消退峰值医疗需求 峰值感染人数 × 重症率6. 模型局限性与改进方向虽然SEIR模型功能强大但也有其局限性。在实际项目中我发现以下几个常见问题需要特别注意数据质量挑战真实疫情数据往往存在报告延迟、检测偏差等问题。曾有一个项目模型预测与官方数据严重不符后来发现是因为检测能力有限导致大量病例未被统计。行为反馈缺失传统模型假设参数恒定但实际上人们看到病例增加后会自发加强防护。可以通过时变参数或行为反馈循环来改进def beta_with_behavioral_feedback(t, I): base_beta 0.4 fear_factor 1 / (1 np.exp(0.1*(I*N - 500))) # 当感染超过500人时行为改变 return base_beta * fear_factor空间异质性大城市中不同区域的接触模式差异很大。一个实用的解决方案是将城市划分为多个区域分别建模后加入区域间迁移项。随机性因素确定性模型无法体现超级传播者等随机事件。可以考虑转向随机微分方程或基于智能体的建模(ABM)# 伪代码基于智能体的简单实现 class Person: def __init__(self): self.state S # S, E, I, or R self.days_infected 0 def interact(self, others, beta): if self.state I: for other in random.sample(others, kpoisson(beta)): if other.state S: other.state E def update(self, sigma, gamma): if self.state E and random.random() sigma: self.state I elif self.state I: self.days_infected 1 if random.random() gamma: self.state R对于希望进一步探索的开发者推荐以下扩展方向加入年龄分层结构不同年龄段接触率和重症率不同考虑无症状感染者和有症状感染者的差异引入疫苗接种状态和疫苗有效性衰减结合网络科学构建接触网络模型