1. 从零开始理解SI模型为什么它适合描述新冠早期传播如果你刚开始接触数学建模或者刚学Python不久看到“SI模型”这个名词可能会有点发怵。别担心这可能是你能接触到的最简单、也最直观的传染病模型之一。我第一次用它来分析问题还是好几年前研究一个校园内信息传播的课题它的简洁和有效让我印象深刻。今天我们就用它来复盘新冠疫情早期的一些传播特征你会发现用几十行Python代码就能模拟出一个复杂疫情的宏观趋势这个过程本身就充满了乐趣和启发性。SI模型全称是Susceptible-Infectious模型中文叫“易感者-感染者模型”。它的核心思想简单到极致把人群分成两类一类是可能被传染的“易感者”S另一类是已经患病并且具有传染能力的“感染者”I。这个模型有一个关键假设感染者一旦被感染就终身具有传染性且不会康复、死亡或获得免疫力。你可能会想这听起来和新冠实际情况不符啊确实新冠有康复者也有重症死亡。但这个简化恰恰是SI模型的用武之地——它非常适合用来刻画一种新病毒在完全无免疫人群中的早期爆发阶段。在疫情最初期医疗认知不足、没有疫苗、没有特效药病毒的传播近乎“野蛮生长”康复和死亡相对于庞大的新增感染数来说影响较小。这时候用SI模型来预测短期内感染人数的激增趋势往往能抓住主要矛盾。用Python来实现它意义重大。对于数学建模新手来说这是一个完美的练手项目涉及的数学是简单的微分方程编程是基础的循环和画图但整个流程——从问题理解、模型建立、到代码实现和结果分析——完整走一遍你对“建模”这件事就会有实实在在的把握感。我们不止是写出代码更要弄懂每一行代码背后的数学含义以及模型结果告诉我们关于疫情传播的哪些秘密。2. SI模型的核心拆解那个驱动疫情的微分方程SI模型的全部奥秘都蕴含在一个微分方程里。我们不需要被“微分方程”吓到它本质上描述的是两类人数量随时间变化的“速度”。假设总人口数 (N) 是固定不变的那么在任何时刻 (t)都有 [ S(t) I(t) N ] 其中(S(t)) 是易感者数量(I(t)) 是感染者数量。模型的关键在于定义感染者数量是如何增加的。SI模型认为新感染者的产生速度即感染者数量 (I) 的增长速度 (\frac{dI}{dt})与易感者数量 (S)和感染者数量 (I)的乘积成正比。为什么是乘积这基于一个合理的接触假设新感染的发生需要易感者和感染者发生有效接触。在一个充分混合的人群中易感者和感染者相遇的机会大致与 (S \times I) 成正比。比例系数记为 (\beta)它被称为感染率或有效接触率它综合反映了病毒的传染力、人群的接触频率等因素。因此我们得到了核心的微分方程 [ \frac{dI}{dt} \beta \cdot \frac{S \cdot I}{N} ]这里有一个细节有时方程写作 (\frac{dI}{dt} \beta S I)有时写作 (\frac{dI}{dt} \beta \frac{S I}{N})。两者的区别在于对 (\beta) 的定义。后者将总人口 (N) 纳入分母这样的 (\beta) 是一个“人均接触率”其含义更清晰且通常是一个介于0和1之间的值更方便理解和赋值。本文我们将采用后一种形式。由于 (S N - I)我们可以把方程改写为只关于 (I(t)) 的方程 [ \frac{dI}{dt} \beta \cdot \frac{(N - I) \cdot I}{N} \beta \cdot i \cdot (1 - i) ] 其中 (i I/N) 是感染者在总人口中的比例。这个方程在数学上被称为Logistic 方程它描述的增长曲线是先慢后快再慢最终趋于饱和的S形曲线。这个方程的解析解是可以求出来的 [ I(t) \frac{N \cdot I_0 \cdot e^{\beta t}}{N I_0 \cdot (e^{\beta t} - 1)} ] 或者用感染比例表示 [ i(t) \frac{i_0 \cdot e^{\beta t}}{1 i_0 \cdot (e^{\beta t} - 1)} ] 其中 (I_0) 或 (i_0) 是初始感染者数量或比例。注意虽然我们有解析解但在实际建模中我们更常用数值方法如欧拉法来求解微分方程。原因有二第一对于更复杂的模型如后续的SIR、SEIR解析解可能不存在或极其复杂第二数值求解的过程能让我们更直观地理解模型随时间步进的动态过程这对初学者理解模型机制至关重要。因此本文将重点介绍数值解法。3. 手把手Python实现从方程到动态曲线理论说得再多不如一行代码。我们使用Python中最常用的科学计算库numpy和绘图库matplotlib来实现SI模型的数值模拟。如果你还没安装可以通过pip install numpy matplotlib快速安装。3.1 参数设置与初始化定义一场“虚拟疫情”首先我们需要定义这场模拟疫情的关键参数。这些参数决定了疫情的走势。import numpy as np import matplotlib.pyplot as plt # 1. 模型参数设置 N 10000 # 总人口假设为一个10000人的封闭社区 I0 1 # 初始感染者人数假设疫情从1个“零号病人”开始 beta 0.3 # 感染率日这是一个关键参数它表示一个感染者平均每天能使多少比例的可疑易感者被感染。 # beta0.3意味着在完全易感的环境中一个感染者每天可能造成约0.3个新感染。 days 100 # 模拟总天数 # 2. 初始化数组 S np.zeros(days) # 易感者数量数组 I np.zeros(days) # 感染者数量数组 S[0] N - I0 # 第0天易感者数量 I[0] I0 # 第0天感染者数量这里beta0.3是一个需要重点理解的参数。它不是一个生物学常数而是一个综合了病毒特性如基本再生数R0和社会行为如接触频率、口罩佩戴的“表现参数”。在疫情早期没有干预措施beta值可能较高随着人们警惕性提高或采取隔离beta值会下降。在SI模型中我们假设它恒定这对应着“无干预”的理想情况。3.2 欧拉法迭代模拟疫情每一天的演进接下来我们用欧拉法来数值求解微分方程。欧拉法的思想很简单用当前时刻的变化率来近似估计下一时刻的值。# 3. 欧拉法迭代求解核心循环 for t in range(0, days-1): # 计算当前时刻t的感染人数变化率 dI/dt # 公式dI/dt beta * (S[t] * I[t]) / N dI_dt beta * (S[t] * I[t]) / N # 更新下一时刻t1的感染者人数和易感者人数 # 欧拉公式I(t1) I(t) dI/dt * delta_t。这里时间步长delta_t1天。 I[t1] I[t] dI_dt # 易感者减少的数量等于感染者增加的数量总人口不变 S[t1] S[t] - dI_dt # 一个重要的安全校验防止计算误差导致人数出现负值或超过总人口 I[t1] np.clip(I[t1], 0, N) S[t1] np.clip(S[t1], 0, N)这个循环是模型的心脏。每一天我们根据当天的S和I计算出当天新增的感染者人数dI_dt然后把它加到累计感染者I中同时从易感者S中减去。np.clip函数是一个实用的技巧确保数值计算中不会因为极小的浮点误差出现负数保证模型的物理意义。3.3 可视化结果解读疫情发展曲线模拟完成数据都在数组里但只有画成图我们才能真正理解发生了什么。# 4. 绘制结果 plt.figure(figsize(10, 6)) time np.arange(days) # 时间轴0到99天 plt.plot(time, S, labelSusceptible (易感者), linewidth2, colorblue) plt.plot(time, I, labelInfectious (感染者), linewidth2, colorred) # 标记一些关键点 # 找到感染者人数超过总人口一半的日期 half_infected_day np.where(I N/2)[0] if len(half_infected_day) 0: plt.axvline(xhalf_infected_day[0], colorgray, linestyle--, alpha0.5) plt.text(half_infected_day[0], N*0.6, fDay {half_infected_day[0]}: 50% Infected, rotation90) plt.xlabel(Time (Days)) plt.ylabel(Number of People) plt.title(SI Model Simulation of COVID-19 Early Spread\n(N{}, I0{}, β{}).format(N, I0, beta)) plt.legend() plt.grid(True, alpha0.3) plt.show()运行这段代码你会得到一张经典的SI模型曲线图。红色感染曲线是一条典型的S形Logistic增长曲线初期缓慢增长期感染者I很少虽然每个感染者传染力强beta固定但能接触到的易感者S有限所以新增缓慢。对应疫情刚出现未被察觉的阶段。中期指数增长期感染者达到一定数量有大量的易感者可供传染每天新增病例数曲线斜率越来越大疫情呈现爆发态势。这对应了疫情在社区中快速扩散的阶段。后期饱和增长期易感者S所剩无几虽然感染者很多但几乎没有新的可传染对象每天新增病例数逐渐减少至零。最终所有人都被感染I趋近于NS趋近于0。这揭示了在无干预、无免疫的SI假设下一种高传染性疾病的必然结局——群体免疫通过极高比例的感染达成而非疫苗。蓝色易感者曲线则是一条从N开始平滑下降至0的曲线。两条曲线在任意时刻的和都等于N。4. 玩转参数beta与I0如何影响疫情走势一个模型的实用性很大程度上取决于我们对参数的理解和把握。在SI模型中beta感染率和I0初始感染者是影响结果最直接的两个参数。我们可以通过简单的参数扫描来直观感受它们的影响力。4.1 感染率beta决定疫情发展的“油门”beta值直接控制了病毒传播的速度。让我们模拟不同beta值下的情况。# 定义一组不同的beta值 beta_values [0.1, 0.2, 0.3, 0.5] plt.figure(figsize(12, 8)) for beta in beta_values: # 重置初始条件 S np.zeros(days) I np.zeros(days) S[0] N - I0 I[0] I0 # 重新模拟 for t in range(0, days-1): dI_dt beta * (S[t] * I[t]) / N I[t1] I[t] dI_dt S[t1] S[t] - dI_dt # 绘制感染者曲线 plt.plot(time, I, labelfβ{beta}, linewidth2) plt.xlabel(Time (Days)) plt.ylabel(Infectious (I)) plt.title(Impact of Infection Rate (β) on Epidemic Spread) plt.legend() plt.grid(True, alpha0.3) plt.axhline(yN, colorblack, linestyle:, alpha0.5, labelTotal Population (N)) plt.show()从图中你可以清晰地看到β0.1疫情发展非常缓慢100天内远未达到全员感染。这说明即使病毒存在如果传播效率很低比如因为人口密度低、卫生习惯好疫情可能只会缓慢蔓延甚至自然消亡在更复杂的模型中。β0.3这是我们基准模拟的情况大约在60-70天时超过一半人口感染。β0.5疫情发展极为迅猛在30天左右就感染了超过一半人口很快达到饱和。实操心得在真实建模中beta是需要通过早期疫情数据来“反演”或“拟合”的关键参数。你可以收集疫情头几周的确诊数据然后调整beta值让你的模型曲线尽可能贴合真实数据。这个过程叫“参数估计”是数学建模连接理论与现实的关键一步。4.2 初始感染者I0引爆疫情的“火星”I0代表了疫情开始的规模。是1个输入病例还是10个、100个这决定了疫情起跑的“起跑线”。I0_values [1, 5, 20, 100] beta_fixed 0.25 plt.figure(figsize(12, 8)) for I0 in I0_values: S np.zeros(days) I np.zeros(days) S[0] N - I0 I[0] I0 for t in range(0, days-1): dI_dt beta_fixed * (S[t] * I[t]) / N I[t1] I[t] dI_dt S[t1] S[t] - dI_dt plt.plot(time, I, labelfI0{I0}, linewidth2) plt.xlabel(Time (Days)) plt.ylabel(Infectious (I)) plt.title(Impact of Initial Infectives (I0) on Epidemic Timeline (β{}).format(beta_fixed)) plt.legend() plt.grid(True, alpha0.3) plt.show()观察结果I01从单个病例开始需要一段较长的“潜伏期”才会进入爆发期。I0100起点很高几乎直接从指数增长期开始疫情曲线陡峭达到饱和的时间大大提前。这个实验说明了早期发现和隔离减小I0对于延缓疫情爆发、为应对争取时间至关重要。即使传播力beta不变控制初始感染规模也能显著改变疫情的时间线。5. 超越基础当SI模型遇到真实世界数据只用模拟数据自娱自乐是不够的。一个合格的建模者总想用模型去触碰真实世界。我们可以尝试用SI模型去拟合新冠疫情早期的真实数据尽管我们知道它有很多局限但这本身是一个极好的学习过程。5.1 数据获取与预处理我们可以从公开数据源如约翰斯·霍普金斯大学CSSE数据集、Our World in Data等获取某个国家或地区疫情早期的累计确诊数据。这里为了演示我们假设手头有一份某城市疫情头30天的简化数据。# 假设我们获取的某城市疫情早期每日累计确诊数据虚构用于演示 # 数据格式[第0天, 第1天, ..., 第29天] real_cumulative_cases np.array([1, 2, 3, 5, 8, 13, 20, 30, 45, 65, 95, 140, 200, 290, 420, 600, 850, 1200, 1700, 2400, 3300, 4500, 6000, 7800, 9800, 12000, 14400, 16800, 19200, 21500]) days_real len(real_cumulative_cases) # 将累计数据转换为每日新增数据更符合SI模型dI/dt的概念 real_new_cases np.diff(real_cumulative_cases) # 计算相邻元素的差值 real_new_cases np.insert(real_new_cases, 0, real_cumulative_cases[0]) # 补回第0天 # 假设该城市总人口N_city为500万 N_city 5_000_000 # 初始感染者I0_city设为第0天的累计确诊数 I0_city real_cumulative_cases[0]5.2 手动拟合与参数估计现在我们的目标是找到一个beta值使得SI模型模拟出的每日新增感染序列与真实的real_new_cases最为接近。这里我们采用一个简单直观的方法——手动调整beta观察拟合效果。# 尝试不同的beta值进行模拟 beta_trials [0.15, 0.22, 0.28, 0.35] plt.figure(figsize(14, 10)) for idx, beta_try in enumerate(beta_trials): # 运行SI模型 S_sim np.zeros(days_real) I_sim np.zeros(days_real) new_cases_sim np.zeros(days_real) # 用于存储模拟的每日新增 S_sim[0] N_city - I0_city I_sim[0] I0_city new_cases_sim[0] 0 # 第一天的新增由后续计算得出 for t in range(0, days_real-1): dI_dt beta_try * (S_sim[t] * I_sim[t]) / N_city I_sim[t1] I_sim[t] dI_dt S_sim[t1] S_sim[t] - dI_dt new_cases_sim[t1] dI_dt # 记录每日新增 # 绘制子图 plt.subplot(2, 2, idx1) plt.plot(range(days_real), real_new_cases, o-, labelReal New Cases (Fictional), linewidth2, markersize4) plt.plot(range(days_real), new_cases_sim, s--, labelfSI Sim (β{beta_try}), linewidth2, markersize4) plt.xlabel(Days) plt.ylabel(New Cases per Day) plt.title(fβ {beta_try}) plt.legend() plt.grid(True, alpha0.3) # 计算一个简单的误差度量均方根误差 (RMSE) rmse np.sqrt(np.mean((real_new_cases - new_cases_sim)**2)) plt.text(0.5* days_real, 0.8*max(real_new_cases), fRMSE: {rmse:.0f}, fontsize10, bboxdict(boxstyleround, facecolorwheat, alpha0.5)) plt.tight_layout() plt.show()通过对比四个子图我们可以发现β0.15模拟曲线增长太慢远低于真实数据RMSE误差很大。β0.28模拟曲线在前期和中期与真实数据贴合得相对较好但后期真实数据增长放缓可能由于干预措施或数据报告延迟而SI模型因无免疫机制仍在快速增长导致后期偏离。β0.22可能是一个折中的拟合值在整个30天内的平均误差较小。注意这是一个非常粗略的手动拟合。在正式研究中我们会使用最小二乘法、极大似然估计等优化算法让计算机自动寻找最优的beta值。但手动尝试的过程能帮你建立对参数敏感性的直觉。5.3 分析拟合结果与模型局限即使我们找到了一个看似不错的beta比如0.22也必须清醒认识到SI模型在此处的严重局限忽略了康复与死亡真实疫情中感染者会康复或死亡从而退出传染链变为Recovered或Removed。这会导致活跃感染者数量I的增长不如SI模型预测的那么快。这也是为什么后期SI模拟值会高于真实新增病例的一个原因。忽略了干预措施在疫情发展过程中政府会采取封控、社交隔离、戴口罩等措施这些都会显著降低有效接触率beta。而我们的模型假设beta恒定。假设人群完全均匀混合模型假设任何易感者和感染者相遇的概率相同忽略了年龄结构、社交网络、空间分布等异质性。所以用SI模型拟合真实数据其目的不在于获得精准的长期预测而在于定性理解早期增长趋势验证疫情初期是否呈现近似指数增长的特征。粗略估计早期传播力通过拟合头一两周的数据可以反推出一个初始的、无干预状态下的beta值这有助于评估病毒的原始传染力。作为教学和更复杂模型的基石理解SI是迈向SIR、SEIR等更现实模型的第一步。6. 从SI到SIR模型的自然演进与Python实现展望当你理解了SI模型自然会想到它的不足并想知道如何改进。最直接、最重要的改进就是加入“康复者”Recovered R群体这就是经典的SIR模型。在SIR模型中人群被分为三类易感者(S)、感染者(I)、康复者(R)。感染者以一定的速率 (\gamma)康复率转变为康复者康复者通常被认为具有免疫力不再被感染。其微分方程组为 [ \begin{align*} \frac{dS}{dt} -\beta \frac{S I}{N} \ \frac{dI}{dt} \beta \frac{S I}{N} - \gamma I \ \frac{dR}{dt} \gamma I \end{align*} ]这个小小的改动增加了一个 (-\gamma I) 项和一个新方程却让模型发生了质变疫情可能不会感染所有人因为康复者获得了免疫他们为易感者提供了一种间接保护。当足够多的人免疫后病毒就找不到足够的易感者来传播疫情就会自然消退。这引出了“群体免疫阈值”的概念。基本再生数R0SIR模型可以自然地推导出 (R_0 \beta / \gamma)。当 (R_0 1) 时疫情会扩散当 (R_0 1) 时疫情会逐渐消失。这成为了衡量传染病传播能力的一个核心指标。用Python实现SIR模型代码结构与我们刚才写的SI模型几乎一样只是多了一个状态变量R和相应的更新逻辑# SIR模型模拟示例伪代码/框架 def simulate_sir(N, I0, beta, gamma, days): S np.zeros(days) I np.zeros(days) R np.zeros(days) S[0] N - I0 I[0] I0 R[0] 0 for t in range(days-1): # 计算变化率 dS_dt -beta * S[t] * I[t] / N dI_dt beta * S[t] * I[t] / N - gamma * I[t] dR_dt gamma * I[t] # 欧拉法更新 S[t1] S[t] dS_dt I[t1] I[t] dI_dt R[t1] R[t] dR_dt # 数值修正 S[t1] np.clip(S[t1], 0, N) I[t1] np.clip(I[t1], 0, N) R[t1] np.clip(R[t1], 0, N) return S, I, R你可以尝试用不同的beta和gamma组合来运行SIR模型观察疫情峰值的出现、最终感染规模等。你会发现通过调整beta代表防控措施强弱和gamma代表医疗水平或病毒毒性模型可以展现出丰富多样的疫情图景这比SI模型更能反映现实世界的复杂性。7. 项目复盘与避坑指南从理论到代码的常见问题走完整个SI模型的构建、实现和分析流程我相信你已经有了不少收获。最后我想分享几个我在教学和实践中学生们最容易踩坑的地方希望能帮你少走弯路。坑一对参数β和γ的物理意义混淆不清问题经常有人把SIR模型中的beta和gamma单位弄混或者不理解gamma1/14是什么意思。解析在SI/SIR模型中时间单位通常是“天”。beta是“感染率”量纲是“每天”表示一个感染者平均每天能传染给多少比例的可疑易感者当S≈N时。gamma是“移出率”在SIR中即康复率量纲也是“每天”。gamma1/14的物理意义是感染者的平均传染期从具有传染性到移出是14天。因为移出过程通常假设为指数衰减平均持续时间就是1/gamma。实操建议在代码注释和报告里务必写明“假设平均传染期为D天则gamma 1/D”。给参数赋予现实意义模型才更有说服力。坑二数值模拟的不稳定性与步长选择问题我们上面用的欧拉法虽然简单但有个缺点如果时间步长我们用的是1天太大或者参数beta很大可能会导致计算出错比如感染人数超过总人口。现象在循环中某一天的dI_dt可能非常大导致I[t1]远大于N虽然我们用np.clip强制修正了但这掩盖了数值方法的不精确。解决方案减小时间步长将1天拆分成更小的步长如0.1天在每个小步长内计算变化。这能显著提高精度。dt 0.1 # 将每天分为10个步长 steps int(days / dt) for i in range(steps-1): dI_dt beta * (S[i] * I[i]) / N * dt # 注意乘以dt I[i1] I[i] dI_dt S[i1] S[i] - dI_dt使用更稳定的数值方法对于兴趣浓厚的同学可以了解并尝试四阶龙格-库塔法RK4它是求解常微分方程的行业标准方法之一精度远高于欧拉法。scipy.integrate库中的odeint或solve_ivp函数内置了更高级的求解器可以直接调用。坑三误用累计数据与新增数据问题在拟合数据时错误地将模型的I(t)累计感染者去拟合真实的“每日新增确诊”数据。解析SI模型输出的I(t)是理论上的累计感染人数。而现实中我们获得的“每日新增确诊”数据更接近模型中的dI/dt每日新增感染人数。此外真实数据还存在报告延迟、检测能力限制等问题并非理论上的“每日新增感染”。实操建议拟合时最好将模型输出的dI/dt序列与真实的“每日新增”序列进行对比。如果只有累计数据可以将其差分后得到近似的新增序列再用于拟合。要时刻意识到模型变量与现实统计指标之间的概念差异。坑四忽视模型的根本假设与适用边界问题拿着SI模型的结果去争论“疫情最终会感染所有人”或者试图用它做长期预测。核心这是最根本的“坑”。所有的模型都是对现实的简化都有其适用范围。SI模型的“无免疫、无康复、无死亡”假设决定了它只适用于描述某种新病原体在完全易感人群中的极早期爆发阶段或者某些信息、谣言在人群中的传播读完信息即“感染”且不会忘记。建议在报告或论文中阐述完模型结果后必须单列一节“模型局限性”Limitations坦诚地说明SI模型的这些简化假设并讨论这些假设在多大程度上偏离了现实以及这种偏离可能如何影响你的结论。这非但不是减分项反而是科学严谨性的体现。通过这个从零开始的SI模型项目你不仅学会了一个经典的传染病模型和它的Python实现更重要的是你体验了数学建模的完整闭环从实际问题出发做出合理假设建立数学模型编写代码求解分析结果并批判性地思考模型的优劣。这个思维框架是应对未来任何建模挑战的宝贵财富。下次当你看到更复杂的SEIR模型、考虑年龄结构的模型、甚至网络传播模型时你会知道它们都是在这个简单而坚固的SI基石上一层层叠加现实世界的复杂性而构建起来的。