凡是参加过数学建模的同学,大概率都有过这样的经历:拿到一道和传染病相关的题目,第一反应是“这个我熟,SIR模型嘛”,于是把方程组一列,用Python跑出两条曲线,写上“预测接下来一个月的感染人数”,觉得大功告成。结果评委给的分数并不高,问题往往不是模型本身错了,而是你根本没有回应题目真正想让你回答的东西——传染病模型在数学建模里从来不是“套公式画图”,它是一套用来支持决策、评估干预、比较策略的定量工具。
这篇文章想跟你聊的,就是传染病模型从“会背公式”到“真正会用”之间那一段路。我会把常见模型的推导逻辑、参数怎么估计、代码怎么写、边界条件在哪、评委真正想看什么,一次讲清楚。内容主要面向准备数学建模竞赛的学生、做数据分析顺便碰到的研究传染病的同学,以及任何想用几个微分方程理解疫情传播规律的人——你不需要很强的数学背景,只要懂一点微积分和Python基础,就能顺着这篇文章把整套流程跑通。
1. 传染病模型在数学建模里的真实角色:不是套公式,而是回答决策问题
1.1 为什么这类题年年有、永远不过时
传染病模型的题目几乎每年都会换着花样出现在各类数学建模竞赛中,因为它天然具备一个优秀建模题的要素:有明确的人群状态划分、有可观测的数据、有政策干预的讨论空间、还有足够的扩展性。你能从最简单的SIR一路做到带隔离、带疫苗、带年龄结构的复杂模型,难度上下限都很大。
但另一个原因是,这类题背后永远站着一个真实的问题:当一个传染病出现时,管理者需要知道现在到底有多严重、接下来会发生什么、该不该封锁、该什么时候解封、疫苗覆盖率要到多少才能形成群体免疫。所有这些问题都能转换成数学语言,而数学建模竞赛的题目,本质上就是在模拟这样一个决策场景。
所以你会发现,历年赛题很少直接说“请建立一个传染病模型”,而是会包装成“请评估某干预措施的效果”“请预测不同接种策略下的感染规模”“请给出一个最优的防控方案”。如果你脑子里只有SIR模型本身,不考虑题目的决策目标,那么模型做得再精细,也拿不到高分。这是我在辅导学生时反复强调的第一件事:先读懂题目在问什么,再决定模型怎么建。
1.2 从“仓室”说起:模型究竟在刻画什么
传染病模型的起点非常简单:把一群人按健康状态分成几个“仓室”(compartment),然后写出各个仓室之间人口流动的方程。比如SIR模型,就是把人分成三类:易感者S(Susceptible)、感染者I(Infectious)、移出者R(Recovered/Removed),分别代表还没得病的人、正在传染别人的人、已经康复并获得免疫或因病退出传播链的人。
这里的核心假设是人群是均匀混合的,也就是说每个人和其他人接触的概率相同,没有空间结构,也没有年龄差异。这个假设在现实中当然不成立,但它的好处是能把问题压缩成几个常微分方程,让我们先抓住传染过程的主要矛盾。
仓室之间的“流动”不是随便画的,它对应着传染病的自然病程:易感者接触到感染者,以一定速度变成感染者;感染者经过一段时间后康复或隔离,变成移出者。这个流动速度由两个关键参数控制——接触率β和恢复率γ。整个模型的精髓就在这两个参数上,后面我会详细展开。
2. SI、SIR、SEIR:四个模型的推导逻辑与选型依据
2.1 SI和SIS:两类最简单的“无免疫”模型
从最简单的情况说起。如果一个人得了病之后不会康复,或者病程短到可以忽略康复,那就只需要考虑S和I两个仓室,这就是SI模型。
假设总人口为N,其中S和I满足S+I=N。每个感染者每天接触一定数量的人,其中易感者所占比例为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/γ共同作用的结果。一个传染病如果R0=3,通常有两种可能:β高但病程短,或者β不算高但病程特别长。两种情况下防控策略完全不一样——前者要减少接触,后者要及时发现隔离。理解了这一点,你就不会在建模时只盯着一个参数了。
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?答案不应该靠感觉,而是看三个问题。第一,题目里是否明确提到了潜伏期或者无症状传播;第二,你手头的数据能否识别出潜伏期的存在;第三,模型的结论是否会对“是否存在潜伏期”这个假设敏感。
我整理了一张选型对照表,你可以在建模时直接参考:
| 模型 | 仓室 | 适用情形 | 关键参数 | 典型结论 |
|---|---|---|---|---|
| SI | S→I | 不康复、短时程的快速传播过程 | β | 最终全部感染 |
| SIS | S→I→S | 无免疫、可反复感染 | β, γ | 地方性流行或清除 |
| SIR | S→I→R | 一次感染终生免疫 | β, γ | 总感染人数与峰值时间 |
| SEIR | S→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_eval=t, method='RK45', 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] # beta=0.3, gamma=0.1, I0取首个数据点 popt, pcov = curve_fit( lambda t, beta, gamma, I0_: fit_new_cases(t, beta, gamma, I0_, S0), t_data, confirmed_new, p0=p0, 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才是报告病例数,k<1。这样可以解释很多“数据不够”的现象。
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.09,R0=3.56,然后呢?你要告诉读者,R0为3.56意味着在没有干预的情况下,疫情处于快速扩散状态,控制它需要把有效接触率至少下降多少,而这可以通过什么措施实现。
第四,也是最要命的,模型结论与题目问题脱节。论文洋洋洒洒写了十几页,却没有一句“根据模型分析,我们建议……”的建议。记住,数学模型是用来回答问题的,不是用来展示数学技巧的。每一章的内容都应该最终指向题目提出的那个决策问题。
我在实际带比赛和评审论文的过程中发现,真正拉开差距的往往不是模型的复杂度,而是对模型假设的清醒认识、对数据的细致处理、以及对结果的合理解释。把基础模型的每一个环节打磨清楚,比堆砌一个看起来高级但没人真正理解的模型要有效得多。传染病模型尤其如此——它看起来门槛低,但几乎所有“高级操作”都是在SIR这个基础上加加减减,地基打不牢,往上是盖不了楼的。
最后给你一个可以直接用的练习路径:找一个真实的历史疫情数据,先用SIR模型拟合,然后分段拟合体现干预影响,接着扩展成SEIR模型,最后做敏感性和多情景分析。把这套流程走一遍,比看十篇优秀论文都有用。等你做完这个练习,再回头看任何传染病建模题目,你都会觉得心里有底。