ARTICLE DETAIL

资讯详情

深耕郑州网站建设与运营推广的一线实战洞察。

从SIR到SEIR:传染病建模的推导、参数估计与Python实战

从SIR到SEIR:传染病建模的推导、参数估计与Python实战 凡是参加过数学建模的同学大概率都有过这样的经历拿到一道和传染病相关的题目第一反应是“这个我熟SIR模型嘛”于是把方程组一列用Python跑出两条曲线写上“预测接下来一个月的感染人数”觉得大功告成。结果评委给的分数并不高问题往往不是模型本身错了而是你根本没有回应题目真正想让你回答的东西——传染病模型在数学建模里从来不是“套公式画图”它是一套用来支持决策、评估干预、比较策略的定量工具。这篇文章想跟你聊的就是传染病模型从“会背公式”到“真正会用”之间那一段路。我会把常见模型的推导逻辑、参数怎么估计、代码怎么写、边界条件在哪、评委真正想看什么一次讲清楚。内容主要面向准备数学建模竞赛的学生、做数据分析顺便碰到的研究传染病的同学以及任何想用几个微分方程理解疫情传播规律的人——你不需要很强的数学背景只要懂一点微积分和Python基础就能顺着这篇文章把整套流程跑通。1. 传染病模型在数学建模里的真实角色不是套公式而是回答决策问题1.1 为什么这类题年年有、永远不过时传染病模型的题目几乎每年都会换着花样出现在各类数学建模竞赛中因为它天然具备一个优秀建模题的要素有明确的人群状态划分、有可观测的数据、有政策干预的讨论空间、还有足够的扩展性。你能从最简单的SIR一路做到带隔离、带疫苗、带年龄结构的复杂模型难度上下限都很大。但另一个原因是这类题背后永远站着一个真实的问题当一个传染病出现时管理者需要知道现在到底有多严重、接下来会发生什么、该不该封锁、该什么时候解封、疫苗覆盖率要到多少才能形成群体免疫。所有这些问题都能转换成数学语言而数学建模竞赛的题目本质上就是在模拟这样一个决策场景。所以你会发现历年赛题很少直接说“请建立一个传染病模型”而是会包装成“请评估某干预措施的效果”“请预测不同接种策略下的感染规模”“请给出一个最优的防控方案”。如果你脑子里只有SIR模型本身不考虑题目的决策目标那么模型做得再精细也拿不到高分。这是我在辅导学生时反复强调的第一件事先读懂题目在问什么再决定模型怎么建。1.2 从“仓室”说起模型究竟在刻画什么传染病模型的起点非常简单把一群人按健康状态分成几个“仓室”compartment然后写出各个仓室之间人口流动的方程。比如SIR模型就是把人分成三类易感者SSusceptible、感染者IInfectious、移出者RRecovered/Removed分别代表还没得病的人、正在传染别人的人、已经康复并获得免疫或因病退出传播链的人。这里的核心假设是人群是均匀混合的也就是说每个人和其他人接触的概率相同没有空间结构也没有年龄差异。这个假设在现实中当然不成立但它的好处是能把问题压缩成几个常微分方程让我们先抓住传染过程的主要矛盾。仓室之间的“流动”不是随便画的它对应着传染病的自然病程易感者接触到感染者以一定速度变成感染者感染者经过一段时间后康复或隔离变成移出者。这个流动速度由两个关键参数控制——接触率β和恢复率γ。整个模型的精髓就在这两个参数上后面我会详细展开。2. SI、SIR、SEIR四个模型的推导逻辑与选型依据2.1 SI和SIS两类最简单的“无免疫”模型从最简单的情况说起。如果一个人得了病之后不会康复或者病程短到可以忽略康复那就只需要考虑S和I两个仓室这就是SI模型。假设总人口为N其中S和I满足SIN。每个感染者每天接触一定数量的人其中易感者所占比例为S/N那么单位时间内新增的感染人数就正比于感染者人数和易感者比例的乘积。写成方程就是[ \frac{dS}{dt} -\beta \frac{S I}{N}, \quad \frac{dI}{dt} \beta \frac{S I}{N} ]这里的β是有效接触率表示一个感染者每天能传染给多少个易感者。SI模型解出来的曲线是经典的S形增长曲线感染人数最终会趋向于总人口N因为没有任何人恢复疫情只会一路蔓延到最后所有人感染。现实中很少有完全对应的场景它可以用来描述某些不产生免疫、或病程极短的急性感染在极短时间内的传播但更多时候是作为教学模型存在。SIS模型则多加了一个恢复项。感染者康复后会回到易感者仓室也就是“得了还能再得”比如普通感冒就接近这个模式。方程变为[ \frac{dS}{dt} -\beta \frac{S I}{N} \gamma I, \quad \frac{dI}{dt} \beta \frac{S I}{N} - \gamma I ]这里γI表示单位时间内康复的人数γ的倒数1/γ就是平均感染期。比如γ0.2意味着平均感染期是5天每天有20%的感染者康复。SIS模型会出现两种不同命运如果β/γ小于某个阈值感染人数会逐渐归零如果超过阈值感染人数会稳定在一个不为零的水平形成地方性流行。这个阈值就是后面要重点说的基本再生数R0的雏形。2.2 SIR模型从微分方程组到再生数SIR模型是竞赛中最常用的模型它在SI的基础上加了一个R仓室感染者康复后进入R且不再被感染。假定总人口N保持不变模型写为[ \frac{dS}{dt} -\beta \frac{S I}{N} ][ \frac{dI}{dt} \beta \frac{S I}{N} - \gamma I ][ \frac{dR}{dt} \gamma I ]很多新手第一次看到这三个方程觉得不过如此但这里有个非常深刻的点从方程中可以直接推导出传染病的“爆发条件”。看dI/dt这一项疫情要扩散需要感染者数量在初期是增加的也就是[ \beta \frac{S}{N} - \gamma 0 ]在疫情刚爆发时几乎所有人都是易感者S/N约等于1于是条件变成β - γ 0也就是β/γ 1。这个无量纲比值就是基本再生数R0它代表在一个完全易感的人群中一个感染者平均能传染给多少人。R0大于1疫情扩散R0小于1疫情自然消退。R0不是“能传染几个人”那种简单说法它是接触率β和病程1/γ共同作用的结果。一个传染病如果R03通常有两种可能β高但病程短或者β不算高但病程特别长。两种情况下防控策略完全不一样——前者要减少接触后者要及时发现隔离。理解了这一点你就不会在建模时只盯着一个参数了。2.3 SEIR模型加入潜伏期后发生了什么SIR模型最大的短板是它默认感染者从“被感染”的那一刻起就具有传染性。但现实中有大量传染病存在潜伏期潜伏期内没有症状、也不传染或传染性很弱。于是SEIR模型在S和I之间插入了一个E仓室Exposed暴露者它代表那些已经感染但尚未具有传染能力的人。[ \frac{dS}{dt} -\beta \frac{S I}{N} ][ \frac{dE}{dt} \beta \frac{S I}{N} - \sigma E ][ \frac{dI}{dt} \sigma E - \gamma I ][ \frac{dR}{dt} \gamma I ]新增的σ是潜伏期转阳率1/σ就是平均潜伏期。模型整体结构仍然是“S→E→I→R”的单向流。潜伏期E这一项的价值在于它让模型的预测曲线相对SIR会更平缓往后推迟而且追踪“有多少人正在潜伏期”对制定隔离策略非常关键——因为潜伏期的人无法通过症状筛查出来这就意味着单靠症状监测是不够的。SEIR还可以继续扩展比如加入无症状感染者、加出生死亡、加隔离仓室、加疫苗接种甚至把人群按照年龄分层。竞赛中到底做到多复杂要看你手头数据能支撑到什么程度。模型不是越复杂越好参数太多而数据太少结果就是过拟合这一点在第5章会专门讲。2.4 到底该用哪个模型一张表和三个判断标准很多同学在此纠结用SIR还是SEIR答案不应该靠感觉而是看三个问题。第一题目里是否明确提到了潜伏期或者无症状传播第二你手头的数据能否识别出潜伏期的存在第三模型的结论是否会对“是否存在潜伏期”这个假设敏感。我整理了一张选型对照表你可以在建模时直接参考模型仓室适用情形关键参数典型结论SIS→I不康复、短时程的快速传播过程β最终全部感染SISS→I→S无免疫、可反复感染β, γ地方性流行或清除SIRS→I→R一次感染终生免疫β, γ总感染人数与峰值时间SEIRS→E→I→R存在潜伏期且潜伏期不传染β, σ, γ潜伏期规模与延迟效应SEIR干预增加隔离/疫苗仓室评估防控策略多参数不同策略下的效果对比判断标准的第四点是“能不能用数据把参数估计出来”。如果题目只给了累计确诊和每日新增你可以识别出γ和β但很难把σ识别得准因为观察数据不直接包含潜伏期信息。这种情况下强行用SEIR反而会让拟合结果极不稳定。不如从SIR入手把基准结论做扎实再在灵敏度分析里说明加入潜伏期会怎样改变结论。这既严谨又稳妥评委挑不出大毛病。3. 让模型真正开口说话参数估计与数据对齐的实操方法3.1 参数β、γ的业务含义与取值范围模型建好了参数从哪来这可能是竞赛中卡住最多人的地方。β和γ不是随便拍脑袋填的它们可以从文献中找到参考值也可以从数据中拟合出来。先说怎么理解这两个参数的量级。γ比较好办它直接对应病程。如果平均感染期是10天那么γ≈0.1/天如果平均感染期是5天γ≈0.2/天。这个数据通常来自医学文献你写论文时可以引用。更麻烦的是β它受病毒本身、人口密度、行为习惯、防控强度等多重因素影响几乎不可能从文献里直接抄一个数。所以实践中普遍的做法是用R0和γ的关系反推β即βR0×γ然后让R0在一个合理范围内做扫描。比如某传染病R0在2到4之间病程10天γ0.1那β就在0.2到0.4之间。你自己跑代码时可以用这个范围作为曲线拟合的初值或约束。这样做的好处是参数有明确意义后续做敏感性分析也方便。3.2 用最小二乘做参数拟合的完整链路当你已经有了每日新增确诊或累计确诊数据最常见的参数估计方法是最小二乘。思路很朴素给定一组β、γ和初始感染人数I0用数值方法解SIR方程得到预测的每日新增或累计值然后和真实数据计算残差平方和不断调节参数让残差最小。具体操作链路分四步。第一步确定目标函数。如果数据是每日新增确诊那对应的模型输出是单位时间内从S流入I的人数也就是βSI/N。如果数据是累计确诊那对应的是N-S(t)也就是已经被感染过的总人数。很多同学算出来的预测值和数据对不上就是在这里搞混了。第二步给出参数的合理初值。初值别乱设用前面说的R0范围推β用病程推γ。curve_fit这类工具虽然是迭代优化但初值差太远很容易收敛到局部最优甚至发散。第三步跑优化。代码可以用scipy.optimize的curve_fit也可以自己写scipy.optimize.minimize。前者方便后者更灵活可以同时对多个参数施加约束。第四步检查拟合效果。不要只看R²有多大还要看残差是不是均匀分布的。如果残差有明显的趋势性比如刚开始拟合得很好后面全部偏离那说明模型结构本身有问题比如忽视了干预措施导致β随时间变化。这种情况下R²再高也不能说明模型可靠。3.3 统计口径差异一个容易被忽视的致命细节参数估计中最容易翻车的其实不是数学而是数据口径。同样是“每日新增”不同渠道可能含义不同是“当日检测阳性人数”还是“当日出现症状的人数”是“本地感染”还是“包含输入病例”是“当日通报”还是“按发病日期回溯”这直接决定了你该用哪个模型输出去拟合。举个例子当日通报的新增病例往往存在周末效应和报告延迟数据序列会出现周期性波动。如果你拿原始通报数据直接拟合得到的参数会有明显的虚假波动。常见的处理手段是取7日移动平均或者用“按发病日期”统计的序列。在做数学建模题时有时题目不会直接给你干净的数据你需要自己在数据预处理阶段把这些细节说清楚并在论文里交代你做了什么处理、为什么这么做。另一个细节是人口基数N。SIR方程里的N会影响传播项βSI/N所以如果你把N取错了数量级拟合出来的β也会跟着错。在竞赛题中研究区域的人口总数通常是给定的如果没给就需要你查资料并明确标注数据来源这也是评委考察信息检索能力的一部分。4. 用Python完整复现一个拟合案例从原始数据到图表4.1 环境准备与数值求解器选型我自己做这类分析时用的是Python的SciPy生态主要是solve_ivp做数值积分curve_fit做参数拟合matplotlib画图。相比自己手写欧拉法或龙格库塔直接用现成求解器更稳定而且自适应步长能避免不少数值问题。初次跑这类代码建议在Jupyter Notebook里做因为需要频繁地调整区间、可视化、检查残差。环境安装没什么特殊要求只要把numpy、scipy、matplotlib装好就行这里不额外展开。用得最多的数值求解器是scipy.integrate.solve_ivp它支持RK45等自适应算法。你不需要懂算法的每一行实现但最好知道一件事用自适应步长方法能够保证在曲线变化剧烈的时候自动缩小步长比用固定步长的欧拉法可靠得多。4.2 核心代码SIR拟合与预测下面是一段可以直接改数据就跑的SIR拟合代码。假设数据是一个numpy数组confirmed_new表示每日新增确诊人数已经按7日移动平均处理过研究区域人口为N。import numpy as np from scipy.integrate import solve_ivp from scipy.optimize import curve_fit import matplotlib.pyplot as plt def sir_ode(t, y, beta, gamma): S, I, R y N S I R dS -beta * S * I / N dI beta * S * I / N - gamma * I dR gamma * I return [dS, dI, dR] def fit_new_cases(t, beta, gamma, I0, S0): # 求解SIR返回每日新增感染人数 beta*S*I/N sol solve_ivp( sir_ode, [t[0], t[-1]], [S0, I0, 0.0], t_evalt, methodRK45, args(beta, gamma) ) S, I, R sol.y new_cases beta * S * I / N return new_cases # 准备数据 t_data np.arange(len(confirmed_new)) N 1000000 # 研究区域总人口 S0 N - confirmed_new[0] I0 confirmed_new[0] p0 [0.3, 0.1, I0] # beta0.3, gamma0.1, I0取首个数据点 popt, pcov curve_fit( lambda t, beta, gamma, I0_: fit_new_cases(t, beta, gamma, I0_, S0), t_data, confirmed_new, p0p0, bounds([0.01, 0.01, 1], [2.0, 1.0, N]) ) beta_fit, gamma_fit, I0_fit popt print(f拟合结果: beta{beta_fit:.4f}, gamma{gamma_fit:.4f}, R0{beta_fit/gamma_fit:.3f})这里有一个细节需要提醒solve_ivp传入的t_eval必须是单调递增的数组而且如果你的数据点特别多跑一次sir模型会稍慢curve_fit的迭代次数也会变多。第一次跑建议先用数据的前半段做拟合得到稳定参数后再用后半段做验证。这样既能展示模型泛化能力又能避免用全量数据拟合后无数据可验证的尴尬。4.3 用图表说话如何展示模型结果代码跑通之后论文里需要三张图缺一不可。第一张是“数据vs模型”的拟合图把真实每日新增确诊和模型预测画在同一个坐标系里让评委一眼看到拟合效果。第二张是S、I、R三条曲线随时间的变化图重点展示感染峰值时间和峰值规模。第三张是参数敏感性分析图通常画R0变化时累计感染人数的变化或者画不同干预强度下的新增曲线对比。我强烈建议在这三张图上面下点功夫因为评委看论文时图的权重非常高。图不是越花哨越好而是要信息清楚、坐标轴标注明确、有图例、有对关键事件的标注。比如你可以在图上标出“干预措施实施日”然后对比该节点前后模型预测的变化这是展示模型决策支持价值最直接的方式。代码层面有两点经验一是matplotlib的中文显示问题提前设置字体否则论文里导出图片会出现方框二是保存图片时用矢量格式PDF或者高分辨率PNG保证印刷和缩放质量。5. 老手也会翻车的五个边界条件5.1 封闭人群假设与人口流动SIR模型默认研究人群是封闭的也就是说没有迁入迁出总人口N恒定。但现实中几乎没有哪个地区是彻底封闭的。在建模竞赛里如果题目给的是某个城市的疫情数据而该城市有大量外来人口那么在疫情初期输入性病例会明显干扰拟合结果。处理方法通常有两种一是把模型的初始条件I0设成大于第一个报告病例数的值用拟合去吸收“存量感染者”的影响二是在模型里显式加入输入项比如在dI/dt上加一个外部输入项Λ(t)只有在题目确实强调输入性风险时才推荐这样做。5.2 β不是常数干预措施如何改写模型这是模型应用层面最大的坑β在现实中根本不是一个常数。戴口罩、保持社交距离、封锁、疫苗全都在改变β。你把一整段疫情数据扔给SIR模型去拟合得到的β只是一个“平均有效接触率”完全没有体现出干预的效果。正确的做法有几种。一是分段拟合按干预措施的实施时间把数据切成几段分别拟合得到不同阶段的β值这样就能量化“封锁使接触率下降了多少”。二是直接让β随时间变化比如设β(t)β0×exp(-kt)用一个衰减函数描述防控不断加强的过程。三是把干预措施作为额外仓室变量显式建模比如增加Q隔离仓室。这些扩展的本质都是承认模型的参数是有业务含义的。评委想看的不是你会不会解SIR方程而是你能不能根据现实背景合理修改模型结构。竞赛论文的加分项往往就体现在这里别人用一个常数β拟合整段数据你把β变成分段函数然后对比分析每一段的下降幅度。5.3 数据延迟与报告误差传染病数据天然存在滞后从感染到出现症状需要几天从出现症状到确诊还需要几天从确诊到通报又需要几天。你手里的“每日新增确诊”并不是“每日真实感染”的同步指标它至少滞后了5到14天。如果不考虑这个滞后模型预测的峰值时间会严重偏离现实。反过来有些同学在拟合时对不上曲线就开始强行调参越调越乱。我见过的比较稳健的做法是在数据预处理阶段明确通报滞后区间然后在模型输出上做一个等长的时间平移来对齐数据。这个方法虽然粗糙但操作简单能有效减少拟合残差并且在论文中说明即可。还有一个容易忽略的问题是漏报。轻症和无症状感染者可能永远不被统计到数据里。这意味着拟合得到的I其实是“被检测到的感染者”不是真实感染人数。如果你发现某段时间新增数据明显偏低可以考虑在模型里加一个检测率参数也就是k×I才是报告病例数k1。这样可以解释很多“数据不够”的现象。5.4 过度拟合与伪预测用微分方程模型做预测和用机器学习模型做预测最大的区别是微分方程模型有强烈的结构约束参数数量很少不容易过度拟合。但如果你开始往模型里加参数——接触率变化、检测率、隔离比例、疫苗生效速度——加到七八个参数以上而数据只有几十个点那模型就开始“记住”数据而不是“理解”数据了。一个非常明显的反面典型是这样拟合优度R²达到0.999看起来完美贴合历史数据但做未来预测时预测曲线要么指数爆炸要么迅速归零。为什么因为参数组合虽然在历史数据上表现很好但在模型结构上完全不稳健微小扰动就会让方程组走向完全不同的状态。怎么避免第一能少加参数就少加。第二做交叉验证用前70%的数据拟合后30%的数据验证。第三报告参数的置信区间。curve_fit返回的pcov就是参数协方差矩阵对角线元素开根号就是标准差。如果你的参数标准差比参数本身还大那基本说明数据信息量不足需要简化模型。这些检验方式在竞赛论文里是非常亮眼的专业细节。5.5 异质性与接触网络最后说一个很多人听过但不知道怎么应对的问题人群不是均匀混合的儿童、成年人、老年人的接触模式差异非常大。同一个城市里有的人一天接触几百人有的人基本不出门。均匀混合的SIR模型相当于假设传染病“平均地”传播这会产生系统性偏差。在竞赛层面上你不需要真的去构建一个个体级别的接触网络但可以做一件事把人群按年龄或活动模式分层建立多组SIR方程组与组之间通过接触矩阵互相感染。这样做之后你会发现模型预测的高峰感染规模、峰值时间都会变而且结论通常更贴近实际。如果题目本身没有要求分层模型往往作为改进方向放在论文的模型扩展部分评委很喜欢看到这种“既有建模基础又有进阶意识”的处理方式。6. 竞赛拿高分的小技巧从模型到报告6.1 敏感性分析与多情景仿真当你已经得到一个拟合好的模型下一步不是急着写结论而是做敏感性分析。简单来说就是系统性改变参数值观察关键输出指标累计感染数、峰值时间、峰值人数如何变化。常用的呈现方式是一个热力图或一组曲线族横轴是β的变化范围纵轴是γ的变化范围颜色代表累计感染人数。多情景仿真则是把敏感性分析包装成决策语言。比如你设定三组情景强干预β下降60%、中等干预β下降30%、弱干预β不变然后分别用模型跑未来的感染曲线对比疫情峰值和结束时间。这种做法的好处是它把“拟合历史数据”升级成了“支撑决策工具”这正是竞赛评审标准里最看重的模型应用价值。6.2 模型验证与残差检查模型验证是很多队伍直接跳过的环节。你要做两件事第一用模型拟合前一段数据预测后一段数据画出预测区间和真实值对比。这是最直观的“模型是否可以外推”的证据。第二检查残差是否像随机噪声。如果你的模型完美拟合了数据但残差呈现明显的周期性——比如每7天一个高峰——那说明模型遗漏了数据本身的周期性特征你需要在预处理中处理周效应而不是改模型。残差检查还有一个额外的好处它能帮你识别异常点。比如某一天的新增病例数突然暴涨那可能是数据口径调整也可能是发生了超级传播事件。在论文里说明这些异常点并解释原因能显著提升报告的可信度。6.3 报告呈现的常见误区最后说几个我在阅卷和辅导中反复看到的报告问题。第一把大段代码贴在正文里。模型建立和求解过程应该用数学公式和文字描述代码放到附录即可。评委关心的是你的建模思路而不是每一行语法。第二图和表没有编号、没有标题、没有在正文中引用。这在数学建模论文中是硬伤直接拉低印象分。第三只报告“拟合出了什么参数”不解释“这些参数意味着什么”。你算出β0.32γ0.09R03.56然后呢你要告诉读者R0为3.56意味着在没有干预的情况下疫情处于快速扩散状态控制它需要把有效接触率至少下降多少而这可以通过什么措施实现。第四也是最要命的模型结论与题目问题脱节。论文洋洋洒洒写了十几页却没有一句“根据模型分析我们建议……”的建议。记住数学模型是用来回答问题的不是用来展示数学技巧的。每一章的内容都应该最终指向题目提出的那个决策问题。我在实际带比赛和评审论文的过程中发现真正拉开差距的往往不是模型的复杂度而是对模型假设的清醒认识、对数据的细致处理、以及对结果的合理解释。把基础模型的每一个环节打磨清楚比堆砌一个看起来高级但没人真正理解的模型要有效得多。传染病模型尤其如此——它看起来门槛低但几乎所有“高级操作”都是在SIR这个基础上加加减减地基打不牢往上是盖不了楼的。最后给你一个可以直接用的练习路径找一个真实的历史疫情数据先用SIR模型拟合然后分段拟合体现干预影响接着扩展成SEIR模型最后做敏感性和多情景分析。把这套流程走一遍比看十篇优秀论文都有用。等你做完这个练习再回头看任何传染病建模题目你都会觉得心里有底。
返回列表