1. 项目概述当新冠疫情遇上SIR模型最近几年大家应该都听过“数学模型预测疫情发展”这类新闻。听起来很高深好像只有专家才能搞懂。其实这背后的核心模型之一——SIR模型其原理并不复杂用我们熟悉的Python就能把它“跑”起来直观地看到疫情是如何传播、达到高峰、然后逐渐平息的。这个项目就是带你用Python亲手搭建一个针对新冠疫情的SIR模型从零开始理解微分方程如何描述人群的动态变化并最终生成一张能讲故事的疫情发展趋势图。无论你是对Python感兴趣的小白还是想参加数学建模比赛找点灵感的同学这个项目都能让你获得“亲手把理论变成可视结果”的成就感。我们将从最基础的模型原理讲起手把手完成代码实现并深入探讨如何调整参数让模型更贴近现实比如隔离措施、疫苗接种这些因素如何影响那条关键的“感染人数曲线”。2. SIR模型核心原理拆解2.1 模型的基本假设与人群划分SIR模型是传染病动力学中最经典的仓室模型之一。它的核心思想非常简单把研究对象的总人口N划分为三个互不重叠的“仓室”。易感者 (Susceptible, S)指那些尚未感染疾病但缺乏免疫力因此有可能被传染的人。在新冠疫情初期几乎所有人都是易感者。感染者 (Infectious, I)指那些已经感染了病毒并且能够将病毒传染给易感者的人。这里模型做了一个简化假设感染者一旦染病就立即具有传染性并且在整个患病期间传染能力不变。移除者 (Removed, R)指那些不会再被感染也不会传染他人的人。这个仓室包含了康复者获得了免疫力和病死者。在SIR模型中他们被同等看待因为从疾病传播动力学的角度看他们都已退出“传播链”。这个划分构成了模型名字的由来S-I-R。模型的关键在于描述这三类人如何随着时间t变化即S(t) I(t) R(t)的动态过程。并且在任何时刻总人口数保持不变S(t) I(t) R(t) N。这是一个非常重要的约束条件我们在编程时需要时刻注意。2.2 微分方程描述变化率的语言模型如何动起来靠的是一组常微分方程。它描述的是每个仓室人数随时间的变化率。我们可以用生活化的类比来理解把S、I、R想象成三个连通的水池水代表人数在水池间流动微分方程就是描述每个水池进水口和出水口水流速度的公式。易感者S的变化率 dS/dt怎么变少只有当易感者接触到感染者时才可能被感染。因此S减少的速度应该与当前易感者的人数S成正比也与当前感染者的人数I成正比接触机会。引入一个比例常数β贝塔称为感染率。所以S的减少速度是 -β * S * I / N。这里除以N是为了将接触概率标准化有时也写作 -β * S * I取决于β的定义是否已包含1/N我们采用更常见的前者。总结dS/dt -β * S * I / N。负号表示S在减少。感染者I的变化率 dI/dt怎么变多易感者被感染后就进入了感染者群体。所以I增加的速度正好就是S减少的速度不考虑因病死亡立即移除的情况即 β * S * I / N。怎么变少感染者不会永远具有传染性。他们要么康复要么死亡都会进入移除者R群体。假设感染者平均持续时间为1/γ天那么每天有一定比例的感染者被移除。这个比例常数γ伽马称为移除率或康复率。I减少的速度是 -γ * I。总结dI/dt β * S * I / N - γ * I。第一项是“进水”新增感染第二项是“出水”移除。移除者R的变化率 dR/dt怎么变多移除者群体的唯一来源就是感染者被移除。所以R增加的速度就是感染者被移除的速度即 γ * I。总结dR/dt γ * I。这一组方程就构成了经典的SIR模型dS/dt -β * S * I / N dI/dt β * S * I / N - γ * I dR/dt γ * I其中β和γ是两个核心参数决定了疫情发展的形态。我们的Python程序本质上就是找一个“计算器”在给定初始值S(0), I(0), R(0)和参数β, γ后把这组方程从第0天算到第100天或任意天数算出每一天的S I R值。注意这里有一个非常重要的概念叫基本再生数R0。它表示在完全易感的人群中一个感染者在其整个传染期内平均能传染的人数。在SIR模型中R0 β / γ。如果R0 1疾病会传播开来如果R0 1疾病会逐渐消失。在新冠疫情初期估算R0对于判断疫情严重性至关重要。3. Python实现从方程到曲线3.1 环境准备与工具选型工欲善其事必先利其器。我们不需要复杂的IDE一个能运行Python的环境加上几个核心库就够了。Python环境确保你安装了Python3.6及以上版本。可以在命令行输入python --version检查。核心库numpy提供高效的数组运算我们的S I R数据将用数组存储。scipy尤其是其中的scipy.integrate.solve_ivp函数它是求解微分方程组的“瑞士军刀”比我们自己写循环更准确、更高效。matplotlib绘图库用来可视化我们的结果。安装这些库非常简单在命令行中使用pip即可pip install numpy scipy matplotlib选择solve_ivp的原因在于SIR模型的微分方程组是“初值问题”给定初始状态求后续变化而solve_ivp实现了多种数值积分方法如龙格-库塔法能稳定、精确地求解这类问题。自己用欧拉法迭代虽然直观但容易累积误差对于学习模型原理后的正式实现推荐直接使用这个工业级工具。3.2 模型定义与参数设定首先我们需要用Python函数的形式来定义我们之前写下的那组微分方程。这个函数将被solve_ivp调用。import numpy as np from scipy.integrate import solve_ivp import matplotlib.pyplot as plt # 1. 定义SIR模型的微分方程组 def sir_model(t, y, beta, gamma, N): 定义SIR模型的微分方程。 参数: t: 时间求解器自动传入 y: 包含[S, I, R]当前值的数组 beta: 感染率 gamma: 移除率 N: 总人口 返回: dydt: [dS/dt, dI/dt, dR/dt]的数组 S, I, R y # 解包当前状态 dS_dt -beta * S * I / N dI_dt beta * S * I / N - gamma * I dR_dt gamma * I return [dS_dt, dI_dt, dR_dt]接下来设定模型参数和初始条件。这些值需要根据实际情况进行估计这里我们使用一组在新冠疫情早期研究中常见的参考值来进行演示。# 2. 设置模型参数和初始条件 N 1000 # 总人口假设为一个1000人的封闭社区 I0 1 # 初始感染者人数假设有1个输入病例 R0 0 # 初始移除者人数 S0 N - I0 - R0 # 初始易感者人数 # 核心参数beta (感染率) 和 gamma (移除率) # 假设平均传染期约为5天则移除率 gamma 1/5 0.2 gamma 0.2 # 基本再生数 R0 是 beta/gamma。假设早期R0约为3.0则 beta R0 * gamma 3.0 * 0.2 0.6 R0_value 3.0 beta R0_value * gamma # 模拟的时间范围 (天) t_start 0 t_end 100 # 模拟100天 t_eval np.linspace(t_start, t_end, 200) # 在100天内均匀取200个时间点用于输出实操心得t_eval参数不是必须的但指定它可以让求解器在我们关心的特定时间点输出结果方便后续绘图。如果不指定求解器会用自己的步长输出的时间点可能不均匀。3.3 方程求解与结果可视化现在万事俱备只差“求解”。我们将调用solve_ivp并把定义好的方程、参数、初始条件都传给它。# 3. 求解微分方程组 # 初始状态向量 initial_state [S0, I0, R0] # 使用solve_ivp求解 # args参数用于向sir_model函数传递额外的参数(beta, gamma, N) solution solve_ivp(sir_model, [t_start, t_end], initial_state, args(beta, gamma, N), t_evalt_eval, methodRK45, # 使用4-5阶龙格-库塔法默认且通常足够精确 dense_outputFalse) # 检查求解是否成功 if not solution.success: print(f求解失败: {solution.message}) else: # 提取结果 t solution.t # 时间点 S solution.y[0] # 易感者人数随时间变化 I solution.y[1] # 感染者人数随时间变化 R solution.y[2] # 移除者人数随时间变化 # 4. 可视化结果 plt.figure(figsize(10, 6)) plt.plot(t, S, label易感者 Susceptible, colorblue, linewidth2) plt.plot(t, I, label感染者 Infectious, colorred, linewidth2) plt.plot(t, R, label移除者 Removed, colorgreen, linewidth2) plt.xlabel(时间 (天), fontsize12) plt.ylabel(人数, fontsize12) plt.title(f新冠疫情SIR模型模拟 (N{N}, $\\beta${beta:.2f}, $\\gamma${gamma:.2f}, $R_0${R0_value:.1f}), fontsize14) plt.legend() plt.grid(True, linestyle--, alpha0.7) plt.tight_layout() plt.show()运行这段代码你将得到一张经典的SIR模型曲线图。红色感染曲线会先上升达到一个峰值后下降最终趋于零。蓝色易感者曲线单调下降绿色移除者曲线单调上升最终易感者和移除者人数之和等于总人口。这张图直观地展示了一场传染病在封闭群体中从发生、爆发到消亡的全过程。4. 深入分析参数影响与模型扩展4.1 关键参数β和γ的敏感性分析模型的“形状”完全由β和γ或者说R0决定。作为建模者我们必须理解改变这些参数意味着什么。感染率β代表了疾病的传播能力与社会接触频率、防控措施强度直接相关。β越大传播越快。模拟实验保持γ0.2不变分别设置β0.3 0.6 0.9对应R01.5 3.0 4.5重新运行模型。你会发现β越大感染峰值越高到来得也越早疫情结束得越快因为易感者迅速耗尽但总感染规模最终R值也越大。移除率γ代表了感染者失去传染性的速度与病程长度或隔离效率相关。γ越大平均传染期越短。模拟实验保持β0.6不变分别设置γ0.1 0.2 0.3对应R06.0 3.0 2.0重新运行。γ越大感染峰值越低疫情曲线越平缓总感染规模越小。这直观地说明了缩短传染期如通过有效隔离和治疗是压平疫情曲线的有效手段。我们可以写一个循环来批量运行并比较不同R0下的情况# 比较不同R0值的影响 gamma_fixed 0.2 R0_values [1.5, 3.0, 4.5] plt.figure(figsize(12, 8)) for R0 in R0_values: beta_temp R0 * gamma_fixed sol solve_ivp(sir_model, [0, 100], [S0, I0, R0], args(beta_temp, gamma_fixed, N), t_evalt_eval) plt.plot(sol.t, sol.y[1], labelf$R_0${R0}, linewidth2) plt.xlabel(时间 (天)) plt.ylabel(感染者人数 I(t)) plt.title(不同基本再生数$R_0$下的感染者曲线对比 ($\\gamma$0.2)) plt.legend() plt.grid(True, linestyle--, alpha0.7) plt.tight_layout() plt.show()4.2 从SIR到更复杂的SEIR模型经典SIR模型假设感染后立即具有传染性这对于像新冠这样有潜伏期的疾病来说是个明显的简化。更贴近现实的模型是SEIR模型它在S和I之间增加了一个潜伏者 (Exposed, E)仓室。潜伏者已感染病毒但尚未出现症状也不具备传染性或传染性很低经过一段平均潜伏期后转为感染者。SEIR模型的微分方程组变为dS/dt -β * S * I / N dE/dt β * S * I / N - σ * E dI/dt σ * E - γ * I dR/dt γ * I其中σ西格玛是潜伏期倒数假设潜伏期平均为1/σ天。用Python实现SEIR模型只需修改模型函数和初始状态向量def seir_model(t, y, beta, sigma, gamma, N): S, E, I, R y dS_dt -beta * S * I / N dE_dt beta * S * I / N - sigma * E dI_dt sigma * E - gamma * I dR_dt gamma * I return [dS_dt, dE_dt, dI_dt, dR_dt] # 参数示例假设平均潜伏期5天则 sigma 1/5 0.2 sigma 0.2 # 初始条件假设有1个潜伏者 E0 1 I0 0 S0 N - E0 - I0 initial_state_seir [S0, E0, I0, 0]加入潜伏期后疫情曲线的上升会更平缓峰值会延迟并可能降低这更符合新冠疫情的观测特征。这个扩展练习能让你深刻理解增加模型细节是如何影响预测结果的。4.3 考虑现实因素隔离、疫苗接种与动态参数要让模型更具现实意义我们可以引入更多机制隔离措施可以简单地将感染率β设置为随时间变化的函数。例如从第30天开始实行严格的社交隔离β值降低。def beta_function(t): if t 30: return 0.6 # 初期传播率 else: return 0.15 # 隔离后传播率 # 在模型方程中将 beta 替换为 beta_function(t)在sir_model函数内部需要接收一个beta_func参数并计算beta beta_func(t)。这能模拟出疫情曲线在干预后迅速被压制的效果。疫苗接种可以视为以一定速率v将易感者S直接转移到移除者R或一个单独的“接种者”仓室。dS/dt -β * S * I / N - v * S dV/dt v * S # (如果增加接种者仓室V)这模拟了通过接种疫苗建立免疫屏障降低易感人群比例从而抑制传播的过程。动态总人口N在长期模拟或病死率较高时可以考虑因病死亡导致总人口缓慢减少。但这会稍微增加模型的复杂度。注意事项每增加一个机制模型就会引入新的参数。参数越多模型“拟合”历史数据的能力可能越强但也更容易“过拟合”并且参数值往往更难从实际数据中准确估计。在数学建模中需要在模型复杂度和实用性之间取得平衡。对于初学者从经典SIR开始理解其核心再逐步增加复杂性是更稳妥的学习路径。5. 实战演练用真实数据校准模型5.1 数据获取与预处理单纯模拟参数意义有限如果我们有一段真实的疫情数据比如某地区每日新增确诊病例数就可以尝试“校准”我们的模型即找到一组β和γ参数使得模型输出的感染曲线与真实数据最吻合。这通常通过参数估计或模型拟合来完成。假设我们从公开数据平台获得了一份简单的每日新增病例数据daily_cases.csv包含date和new_cases两列。import pandas as pd # 读取数据 data pd.read_csv(daily_cases.csv) data[date] pd.to_datetime(data[date]) data data.sort_values(date).reset_index(dropTrue) # 计算累计病例作为移除者R的近似简化假设所有病例最终都移除 # 注意真实数据中R应包括康复和死亡且存在报告延迟这里仅为演示。 data[cumulative_cases] data[new_cases].cumsum() # 确定总人口N和初始条件 N 10_000_000 # 假设该地区人口1000万 I0 data.loc[0, new_cases] # 假设第一天的新增即为初始感染者粗略 R0 0 S0 N - I0 - R05.2 定义损失函数与参数优化我们的目标是调整β和γ让模型预测的每日新增感染人数即dI/dt dR/dt的某种组合严格来说是dI/dt γ*I但常简化为dI/dt的离散积分与真实的每日新增病例数据之间的差距最小。这个差距可以用误差的平方和来衡量。from scipy.optimize import minimize def model_output(params, data, N): 给定参数运行SIR模型并返回模拟的每日新增感染者数 beta, gamma params initial_state [S0, I0, R0] t_span [0, len(data)-1] t_eval np.arange(len(data)) sol solve_ivp(sir_model, t_span, initial_state, args(beta, gamma, N), t_evalt_eval, methodRK45) if not sol.success: return np.inf * np.ones(len(data)) S, I, R sol.y # 计算模拟的每日新增感染者的变化量 移除者的变化量 (近似为新增报告) # 更简单的做法计算I的差分后一天减前一天但这不稳定。常用的是计算每日从S到I的转移人数。 # 这里采用一种近似每日新增 ≈ β * S * I / N 即新感染人数 simulated_new_infections beta * S * I / N return simulated_new_infections def loss_function(params, data, N): 损失函数计算模拟新增与真实新增的均方误差 real_new_cases data[new_cases].values simulated model_output(params, data, N) # 确保模拟结果有效且长度匹配 if np.any(np.isinf(simulated)) or len(simulated) ! len(real_new_cases): return 1e10 # 返回一个很大的数作为惩罚 mse np.mean((simulated - real_new_cases) ** 2) return mse # 设置参数初始猜测值和边界 initial_guess [0.5, 0.1] # [beta, gamma]的初始猜测 bounds [(0.001, 1.0), (0.01, 0.5)] # 给参数设定合理的物理边界 # 执行优化 result minimize(loss_function, initial_guess, args(data, N), boundsbounds, methodL-BFGS-B) if result.success: beta_estimated, gamma_estimated result.x R0_estimated beta_estimated / gamma_estimated print(f优化得到的参数: beta {beta_estimated:.4f}, gamma {gamma_estimated:.4f}) print(f估算的基本再生数 R0 {R0_estimated:.2f}) else: print(参数优化失败:, result.message)5.3 结果对比与模型评估得到最优参数后我们用这组参数重新运行一次模型并将模拟曲线与真实数据画在一起对比。# 使用估计的参数进行最终模拟 final_params (beta_estimated, gamma_estimated) simulated_new model_output(final_params, data, N) # 绘图对比 plt.figure(figsize(12, 6)) plt.plot(data.index, data[new_cases], o-, label真实每日新增病例, markersize4, linewidth1.5) plt.plot(data.index, simulated_new, r--, labelf模型拟合 (β{beta_estimated:.3f}, γ{gamma_estimated:.3f}), linewidth2) plt.xlabel(时间 (天)) plt.ylabel(人数) plt.title(SIR模型拟合新冠疫情每日新增病例数据对比) plt.legend() plt.grid(True, linestyle--, alpha0.7) plt.tight_layout() plt.show() # 也可以绘制S, I, R三条曲线的完整模拟结果 sol_final solve_ivp(sir_model, [0, len(data)-1], [S0, I0, R0], args(beta_estimated, gamma_estimated, N), t_evalnp.linspace(0, len(data)-1, 300)) t_fine sol_final.t S_final, I_final, R_final sol_final.y plt.figure(figsize(10, 6)) plt.plot(t_fine, S_final, b-, label易感者 S) plt.plot(t_fine, I_final, r-, label感染者 I) plt.plot(t_fine, R_final, g-, label移除者 R) plt.xlabel(时间 (天)) plt.ylabel(人数) plt.title(基于拟合参数的SIR模型模拟曲线) plt.legend() plt.grid(True) plt.tight_layout() plt.show()通过对比图你可以直观地看到模型在多大程度上能还原真实的疫情趋势。拟合不可能完美因为真实世界的影响因素远比SIR模型复杂如检测能力波动、防控政策变化、人口流动等但一个好的拟合能帮助我们理解疫情传播的总体动力学特征并对关键参数如R0做出估计。常见问题与排查技巧拟合结果很差曲线完全对不上检查初始条件I0的设置非常关键。如果数据初期病例数很少设为1可能合适如果数据已经处于爆发期可能需要将初始感染者设为更早时间点的累计病例估计值。检查损失函数确保模拟输出与真实数据的量纲和意义对齐。本例中用β * S * I / N近似每日新增这只在特定假设下成立。更严谨的做法是将模型输出的累计感染人数与真实累计数据对比或者使用更复杂的似然函数。尝试不同优化方法/初始值minimize可能陷入局部最优。尝试多组不同的initial_guess或者使用全局优化算法如basinhopping。求解器警告或失败可能是参数值导致方程数值不稳定如β极大。为参数设置合理的bounds约束很重要。尝试调整solve_ivp的容差参数rtol和atol例如设为rtol1e-6, atol1e-8提高计算精度。模型曲线与真实数据存在系统性偏差这可能暗示模型结构本身不足以描述数据。例如真实数据呈现多波峰而SIR模型通常只产生单峰。这时就需要考虑引入时变参数如随时间衰减的β模拟防控加强或更复杂的模型结构如SEIR或考虑无症状感染者的SAIR模型。这个从理论到实践从模拟到拟合的完整流程是数学建模解决实际问题的核心缩影。通过Python我们不仅实现了模型还完成了关键的一步——用数据去验证和修正它。尽管SIR模型是对现实的极大简化但它为我们提供了分析传染病传播的一个强大而清晰的框架。当你下次再看到关于疫情预测的新闻时希望你能想起这几行代码和它们背后的思想那便是这个项目最大的价值。