简介这份文档面向学习传染病动力学、数学建模或微分方程应用的高校学生与研究者系统讲解SIR模型的建模思路与求解方法。内容从Kermack与McKendrick的经典假设出发梳理易感者、感染者、康复者三类人群的划分给出ds/dt、di/dt、dr/dt微分方程组并结合λ1、μ0.3、s(0)0.98、i(0)0.02的算例用MATLAB的ode45进行数值积分展示感染人数在t≈7达到峰值后衰减、易感者比例趋于稳定非零值的规律还涉及相轨线分析与模型局限性讨论。资源包共1个doc文件约360KB以文档形式完整呈现模型假设、构成、数值计算与意义评估便于直接阅读与引用。目前已有46人学习适合作为课程作业、论文建模或公共卫生防控策略分析的参考起点。1. 传染病问题中的SIR模型从微分方程到可运行代码的落地路径很多人第一次接触传染病建模是在一份名为“传染病问题中的SIR模型.doc”的文档里看到三个耦合微分方程然后被要求用欧拉法或者龙格-库塔法数值求解。但真正做过一轮完整落地的人会告诉你把方程敲进代码只是开始参数标定、初值选取、结果验证才是决定模型能不能用的关键。SIR模型把人群分成易感者Susceptible、感染者Infectious、移除者Removed三类用三个常微分方程描述状态转移结构极简却是后续SEIR、SEIRS等扩展模型的地基。这篇文章面向需要用SIR模型做实际预测或教学复现的工程师和数据分析人员从方程推导讲到Python实现、参数估计、避坑排查最后落到一个能直接跑的完整流程。如果你手里正好有类似“传染病问题中的SIR模型.doc”这样的材料想把它变成可运行、可验证的代码下面的内容可以照着走一遍。2. SIR模型的数学骨架与离散化选择2.1 三个方程到底在描述什么SIR模型的核心假设是人群总数N在短期疫情中保持不变忽略出生和自然死亡只考虑感染和恢复两个过程。设S(t)、I(t)、R(t)分别为t时刻三类人群的数量或比例标准形式如下dS/dt -β * S * I / N dI/dt β * S * I / N - γ * I dR/dt γ * Iβ是传染率系数表示单位时间内一个感染者与易感者接触并成功传播的概率乘以接触数γ是移除率等于1/感染期。这里有一个容易被忽略的细节β和γ的量纲必须匹配如果时间单位是天γ就是每天恢复的比例感染期就是1/γ天。基本再生数R0 β/γ这是判断疫情是否扩散的阈值——R0大于1时感染人数会先升后降小于1时单调衰减。很多文档只给方程不给量纲说明导致后面参数估计时单位混乱。我一般会在代码里显式标注每个参数的单位比如β单位是“1/天”γ单位是“1/天”这样后续调参不会翻车。2.2 数值解法欧拉法够不够用SIR模型没有解析解必须数值求解。常见做法是欧拉法、改进欧拉法或四阶龙格-库塔法RK4。欧拉法最简单但步长稍大就会累积误差甚至出现S(t)变成负数这种物理上不可能的结果。RK4精度高但代码稍复杂。对于大多数疫情数据天级别粒度我一般用RK4步长取0.1天再降采样到天兼顾精度和速度。下面是一个用Python实现RK4求解SIR的最小代码块import numpy as np def sir_deriv(t, y, beta, gamma, N): S, I, R y dS -beta * S * I / N dI beta * S * I / N - gamma * I dR gamma * I return np.array([dS, dI, dR]) def sir_rk4(beta, gamma, N, I0, R0, days, dt0.1): # 初值S0 N - I0 - R0 y np.array([N - I0 - R0, I0, R0], dtypefloat) t 0.0 results [(t, *y)] steps int(days / dt) for _ in range(steps): k1 sir_deriv(t, y, beta, gamma, N) k2 sir_deriv(t dt/2, y dt*k1/2, beta, gamma, N) k3 sir_deriv(t dt/2, y dt*k2/2, beta, gamma, N) k4 sir_deriv(t dt, y dt*k3, beta, gamma, N) y y dt * (k1 2*k2 2*k3 k4) / 6 t dt results.append((t, *y)) return np.array(results)逻辑说明sir_deriv返回三个微分方程的右端项sir_rk4按RK4公式迭代。参数dt是内部步长days是模拟总天数I0和R0是初始感染者和移除者数量。注意R0这个变量名和基本再生数R0重名了代码里用R0表示初始移除人数实际项目中建议改成R_init避免混淆。运行后results的每一行是(t, S, I, R)可以直接用matplotlib画曲线。2.3 参数β和γ的估计最小二乘还是网格搜索有了求解器下一步是让模型贴合真实数据。β和γ未知需要从累计确诊或现存感染者数据中反推。常见做法有两种一是用scipy.optimize.minimize做最小二乘拟合二是网格搜索加手动调参。最小二乘更快但目标函数非凸时容易陷入局部最优网格搜索慢但能看清参数空间的全貌。我一般先用网格搜索粗定位再用最小二乘精修。目标函数用预测的I(t)和实际现存感染者数的残差平方和。注意实际数据里“现存感染者”对应I(t)“累计移除”对应R(t)累计确诊对应I(t)R(t)别搞混。from scipy.optimize import minimize import numpy as np def loss(params, data_I, N, I0, R0, days): beta, gamma params if beta 0 or gamma 0: return 1e10 sol sir_rk4(beta, gamma, N, I0, R0, days, dt0.1) # 降采样到天 daily sol[::10] pred_I daily[:, 2] return np.sum((pred_I - data_I) ** 2) # 假设data_I是每天现存感染者数组长度days1 res minimize(loss, x0[0.3, 0.1], args(data_I, N, I0, R0, days), methodNelder-Mead, options{xatol: 1e-6, fatol: 1e-6}) beta_hat, gamma_hat res.x参数说明x0是β和γ的初始猜测一般β取0.2~0.5γ取0.05~0.2。methodNelder-Mead不需要梯度适合这种黑箱目标函数。如果数据噪声大可以加正则项或者改用L-BFGS-B并手动计算梯度。拟合完成后务必检查R0 beta_hat / gamma_hat是否合理如果R0小于1但数据明显在增长说明拟合有问题。3. 用真实数据跑通SIR数据预处理与拟合流程3.1 数据从哪来、怎么清洗做SIR建模数据质量比模型本身更重要。常见数据源包括公开的疫情时间序列、医院上报的每日新增和现存病例。拿到数据后先做三件事对齐时间轴、处理缺失值、把累计量转成现存量。累计确诊减去累计死亡和累计治愈就是现存感染者但很多数据集里治愈和死亡有滞后上报直接相减会得到负数。我一般用滑动平均或者样条插值平滑后再相减。另外SIR模型假设移除者不再感染如果数据里出现“复阳”或者二次感染这个假设就不成立需要考虑SEIRS。对于短期预测比如一个月内SIR通常够用但要在报告里注明假设。3.2 完整拟合脚本从CSV到参数输出下面是一个从CSV文件读取数据、拟合β和γ、输出R0和预测曲线的完整脚本。假设CSV有两列date和active_cases。import pandas as pd import numpy as np from scipy.optimize import minimize import matplotlib.pyplot as plt # 读取数据 df pd.read_csv(cases.csv, parse_dates[date]) df df.sort_values(date).reset_index(dropTrue) data_I df[active_cases].values.astype(float) # 参数设置 N 1_000_000 # 总人口根据实际地区调整 I0 data_I[0] # 初始感染者 R0_init 0 # 初始移除者假设为0 days len(data_I) - 1 # 定义损失函数同上 def loss(params, data_I, N, I0, R_init, days): beta, gamma params if beta 0 or gamma 0: return 1e10 sol sir_rk4(beta, gamma, N, I0, R_init, days, dt0.1) daily sol[::10] pred_I daily[:, 2] return np.sum((pred_I - data_I) ** 2) # 拟合 res minimize(loss, x0[0.3, 0.1], args(data_I, N, I0, R0_init, days), methodNelder-Mead, options{xatol: 1e-6, fatol: 1e-6}) beta_hat, gamma_hat res.x R0 beta_hat / gamma_hat print(fbeta{beta_hat:.4f}, gamma{gamma_hat:.4f}, R0{R0:.2f}) # 预测并画图 sol sir_rk4(beta_hat, gamma_hat, N, I0, R0_init, days, dt0.1) daily sol[::10] plt.plot(df[date], data_I, o, label实际) plt.plot(df[date], daily[:, 2], -, labelSIR预测) plt.legend() plt.show()逻辑说明先读CSV并排序提取现存感染者列。N需要根据实际地区人口设置如果数据是比例而非绝对人数N取1即可。拟合用Nelder-Mead初始值给0.3和0.1。输出R0后画图对比。如果预测曲线和实际偏差大先检查数据是否有异常点再检查N和I0是否合理。参数调整建议如果拟合出的γ对应的感染期明显偏离常识比如超过30天可能是数据里混入了长期阳性病例需要截断或者改用更复杂的模型。β的估计对早期数据敏感可以只用前70%的数据拟合后30%做验证。3.3 模型验证留出法和残差检查拟合完不能只看R0还要做验证。最简单的是留出法用前80%的数据拟合预测后20%计算预测误差。如果误差在可接受范围内比如MAPE小于20%模型可用。另外画残差图如果残差随时间呈现系统性偏差说明模型结构有问题比如需要加入时变β或者干预措施。我一般还会检查S(t)I(t)R(t)是否始终等于N如果数值解出现偏离说明步长太大或者RK4实现有误。这个检查能快速定位低级错误。4. 避坑与排查SIR建模中最容易翻车的五个地方4.1 现象S(t)变成负数曲线出现震荡原因欧拉法步长过大或者RK4实现时dt没有同步更新。SIR方程在I接近0时刚性较强显式方法容易不稳定。解决减小dt到0.01以下或者改用自适应步长的scipy.integrate.solve_ivp。如果必须用固定步长确保RK4的四个斜率都用同一个dt计算。4.2 现象拟合出的R0高得离谱比如大于10原因数据里的现存感染者被低估或者N设置得太小。比如实际人口100万代码里N写成1万β就会被放大。解决核对N的单位和数据单位是否一致。如果数据是每万人中的病例数N就取1。另外检查I0是否用了第一天的数据如果第一天数据是0I0要手动设一个大于0的值。4.3 现象预测曲线比实际晚很多天到达峰值原因γ估计偏小感染期被拉长。常见于数据里包含大量轻症或无症状感染者他们实际移除更快但被计入了长期感染者。解决用累计确诊和累计移除分别拟合或者改用SEIR把潜伏期单独建模。如果坚持用SIR可以对γ加一个下界约束比如gamma 0.05。4.4 现象不同时间段拟合出的β差异巨大原因真实疫情中干预措施戴口罩、限流会改变接触率β不是常数。SIR的常数β假设在长周期内不成立。解决把时间分段每段单独拟合β或者用时变β(t)的扩展模型。如果只是做短期预测用最近两周的数据拟合即可。4.5 现象代码跑得通但结果和文档里的对不上原因文档里的方程可能用了不同的符号约定比如β定义为单位时间接触数而不是传染概率或者R(t)表示累计确诊而不是移除者。解决先用手算一个简单例子验证。取N1000I01β0.3γ0.1手动算前两步再和代码输出对比。这个步骤能省掉后面几小时的排查时间。5. 进阶技巧用SIR做情景分析和干预模拟5.1 把β拆成接触率和传染概率基础SIR里β是一个黑箱参数但实际做干预评估时需要知道改变哪个环节有效。常见做法是把β拆成c人均接触数和p每次接触的传染概率β c * p。这样戴口罩降低p限流降低c可以分别模拟。代码上只需要在sir_deriv里把beta替换成c*p然后对c和p分别扫描。def sir_deriv_intervention(t, y, c, p, gamma, N): S, I, R y beta c * p dS -beta * S * I / N dI beta * S * I / N - gamma * I dR gamma * I return np.array([dS, dI, dR])参数说明c的单位是“接触数/天”p是无量纲概率。如果原来β0.3可以设c10p0.03。干预模拟时把c降到5看峰值如何变化。5.2 用蒙特卡洛评估参数不确定性拟合出的β和γ有置信区间直接拿点估计做预测会低估风险。我一般用蒙特卡洛从拟合的协方差矩阵里采样1000组参数每组跑一次SIR取I(t)的5%和95%分位数作为预测区间。这样报告里能给出“峰值在X到Y之间”的结论比单条曲线可信得多。import numpy as np # 假设res是拟合结果res.hess_inv是逆Hessian近似 params_samples np.random.multivariate_normal(res.x, res.hess_inv, size1000) peak_I [] for beta_s, gamma_s in params_samples: if beta_s 0 or gamma_s 0: continue sol sir_rk4(beta_s, gamma_s, N, I0, R0_init, days, dt0.1) peak_I.append(sol[:, 2].max()) peak_I np.array(peak_I) print(f峰值感染者中位数{np.median(peak_I):.0f}) print(f90%置信区间[{np.percentile(peak_I, 5):.0f}, {np.percentile(peak_I, 95):.0f}])注意res.hess_inv在Nelder-Mead下可能不可用可以改用L-BFGS-B或者用bootstrap重采样。这个技巧在写报告时特别有用能避免“模型预测峰值是1000”这种过于绝对的表述。5.3 一个我常犯的错误早期我做SIR拟合时总想把所有数据点都拟合得很准结果β和γ被少数异常点带偏。后来养成习惯先画散点图肉眼剔除明显离群的天数再拟合。另外SIR模型对早期数据敏感如果第一天只有1个病例I01拟合出的β会偏大。我一般把前3天的数据去掉从第4天开始拟合用前3天估算I0。这个习惯帮我省了很多次返工。希望帮到你。本文还有配套的精品资源点击获取
传染病SIR模型实战：从微分方程到Python代码落地
简介这份文档面向学习传染病动力学、数学建模或微分方程应用的高校学生与研究者系统讲解SIR模型的建模思路与求解方法。内容从Kermack与McKendrick的经典假设出发梳理易感者、感染者、康复者三类人群的划分给出ds/dt、di/dt、dr/dt微分方程组并结合λ1、μ0.3、s(0)0.98、i(0)0.02的算例用MATLAB的ode45进行数值积分展示感染人数在t≈7达到峰值后衰减、易感者比例趋于稳定非零值的规律还涉及相轨线分析与模型局限性讨论。资源包共1个doc文件约360KB以文档形式完整呈现模型假设、构成、数值计算与意义评估便于直接阅读与引用。目前已有46人学习适合作为课程作业、论文建模或公共卫生防控策略分析的参考起点。1. 传染病问题中的SIR模型从微分方程到可运行代码的落地路径很多人第一次接触传染病建模是在一份名为“传染病问题中的SIR模型.doc”的文档里看到三个耦合微分方程然后被要求用欧拉法或者龙格-库塔法数值求解。但真正做过一轮完整落地的人会告诉你把方程敲进代码只是开始参数标定、初值选取、结果验证才是决定模型能不能用的关键。SIR模型把人群分成易感者Susceptible、感染者Infectious、移除者Removed三类用三个常微分方程描述状态转移结构极简却是后续SEIR、SEIRS等扩展模型的地基。这篇文章面向需要用SIR模型做实际预测或教学复现的工程师和数据分析人员从方程推导讲到Python实现、参数估计、避坑排查最后落到一个能直接跑的完整流程。如果你手里正好有类似“传染病问题中的SIR模型.doc”这样的材料想把它变成可运行、可验证的代码下面的内容可以照着走一遍。2. SIR模型的数学骨架与离散化选择2.1 三个方程到底在描述什么SIR模型的核心假设是人群总数N在短期疫情中保持不变忽略出生和自然死亡只考虑感染和恢复两个过程。设S(t)、I(t)、R(t)分别为t时刻三类人群的数量或比例标准形式如下dS/dt -β * S * I / N dI/dt β * S * I / N - γ * I dR/dt γ * Iβ是传染率系数表示单位时间内一个感染者与易感者接触并成功传播的概率乘以接触数γ是移除率等于1/感染期。这里有一个容易被忽略的细节β和γ的量纲必须匹配如果时间单位是天γ就是每天恢复的比例感染期就是1/γ天。基本再生数R0 β/γ这是判断疫情是否扩散的阈值——R0大于1时感染人数会先升后降小于1时单调衰减。很多文档只给方程不给量纲说明导致后面参数估计时单位混乱。我一般会在代码里显式标注每个参数的单位比如β单位是“1/天”γ单位是“1/天”这样后续调参不会翻车。2.2 数值解法欧拉法够不够用SIR模型没有解析解必须数值求解。常见做法是欧拉法、改进欧拉法或四阶龙格-库塔法RK4。欧拉法最简单但步长稍大就会累积误差甚至出现S(t)变成负数这种物理上不可能的结果。RK4精度高但代码稍复杂。对于大多数疫情数据天级别粒度我一般用RK4步长取0.1天再降采样到天兼顾精度和速度。下面是一个用Python实现RK4求解SIR的最小代码块import numpy as np def sir_deriv(t, y, beta, gamma, N): S, I, R y dS -beta * S * I / N dI beta * S * I / N - gamma * I dR gamma * I return np.array([dS, dI, dR]) def sir_rk4(beta, gamma, N, I0, R0, days, dt0.1): # 初值S0 N - I0 - R0 y np.array([N - I0 - R0, I0, R0], dtypefloat) t 0.0 results [(t, *y)] steps int(days / dt) for _ in range(steps): k1 sir_deriv(t, y, beta, gamma, N) k2 sir_deriv(t dt/2, y dt*k1/2, beta, gamma, N) k3 sir_deriv(t dt/2, y dt*k2/2, beta, gamma, N) k4 sir_deriv(t dt, y dt*k3, beta, gamma, N) y y dt * (k1 2*k2 2*k3 k4) / 6 t dt results.append((t, *y)) return np.array(results)逻辑说明sir_deriv返回三个微分方程的右端项sir_rk4按RK4公式迭代。参数dt是内部步长days是模拟总天数I0和R0是初始感染者和移除者数量。注意R0这个变量名和基本再生数R0重名了代码里用R0表示初始移除人数实际项目中建议改成R_init避免混淆。运行后results的每一行是(t, S, I, R)可以直接用matplotlib画曲线。2.3 参数β和γ的估计最小二乘还是网格搜索有了求解器下一步是让模型贴合真实数据。β和γ未知需要从累计确诊或现存感染者数据中反推。常见做法有两种一是用scipy.optimize.minimize做最小二乘拟合二是网格搜索加手动调参。最小二乘更快但目标函数非凸时容易陷入局部最优网格搜索慢但能看清参数空间的全貌。我一般先用网格搜索粗定位再用最小二乘精修。目标函数用预测的I(t)和实际现存感染者数的残差平方和。注意实际数据里“现存感染者”对应I(t)“累计移除”对应R(t)累计确诊对应I(t)R(t)别搞混。from scipy.optimize import minimize import numpy as np def loss(params, data_I, N, I0, R0, days): beta, gamma params if beta 0 or gamma 0: return 1e10 sol sir_rk4(beta, gamma, N, I0, R0, days, dt0.1) # 降采样到天 daily sol[::10] pred_I daily[:, 2] return np.sum((pred_I - data_I) ** 2) # 假设data_I是每天现存感染者数组长度days1 res minimize(loss, x0[0.3, 0.1], args(data_I, N, I0, R0, days), methodNelder-Mead, options{xatol: 1e-6, fatol: 1e-6}) beta_hat, gamma_hat res.x参数说明x0是β和γ的初始猜测一般β取0.2~0.5γ取0.05~0.2。methodNelder-Mead不需要梯度适合这种黑箱目标函数。如果数据噪声大可以加正则项或者改用L-BFGS-B并手动计算梯度。拟合完成后务必检查R0 beta_hat / gamma_hat是否合理如果R0小于1但数据明显在增长说明拟合有问题。3. 用真实数据跑通SIR数据预处理与拟合流程3.1 数据从哪来、怎么清洗做SIR建模数据质量比模型本身更重要。常见数据源包括公开的疫情时间序列、医院上报的每日新增和现存病例。拿到数据后先做三件事对齐时间轴、处理缺失值、把累计量转成现存量。累计确诊减去累计死亡和累计治愈就是现存感染者但很多数据集里治愈和死亡有滞后上报直接相减会得到负数。我一般用滑动平均或者样条插值平滑后再相减。另外SIR模型假设移除者不再感染如果数据里出现“复阳”或者二次感染这个假设就不成立需要考虑SEIRS。对于短期预测比如一个月内SIR通常够用但要在报告里注明假设。3.2 完整拟合脚本从CSV到参数输出下面是一个从CSV文件读取数据、拟合β和γ、输出R0和预测曲线的完整脚本。假设CSV有两列date和active_cases。import pandas as pd import numpy as np from scipy.optimize import minimize import matplotlib.pyplot as plt # 读取数据 df pd.read_csv(cases.csv, parse_dates[date]) df df.sort_values(date).reset_index(dropTrue) data_I df[active_cases].values.astype(float) # 参数设置 N 1_000_000 # 总人口根据实际地区调整 I0 data_I[0] # 初始感染者 R0_init 0 # 初始移除者假设为0 days len(data_I) - 1 # 定义损失函数同上 def loss(params, data_I, N, I0, R_init, days): beta, gamma params if beta 0 or gamma 0: return 1e10 sol sir_rk4(beta, gamma, N, I0, R_init, days, dt0.1) daily sol[::10] pred_I daily[:, 2] return np.sum((pred_I - data_I) ** 2) # 拟合 res minimize(loss, x0[0.3, 0.1], args(data_I, N, I0, R0_init, days), methodNelder-Mead, options{xatol: 1e-6, fatol: 1e-6}) beta_hat, gamma_hat res.x R0 beta_hat / gamma_hat print(fbeta{beta_hat:.4f}, gamma{gamma_hat:.4f}, R0{R0:.2f}) # 预测并画图 sol sir_rk4(beta_hat, gamma_hat, N, I0, R0_init, days, dt0.1) daily sol[::10] plt.plot(df[date], data_I, o, label实际) plt.plot(df[date], daily[:, 2], -, labelSIR预测) plt.legend() plt.show()逻辑说明先读CSV并排序提取现存感染者列。N需要根据实际地区人口设置如果数据是比例而非绝对人数N取1即可。拟合用Nelder-Mead初始值给0.3和0.1。输出R0后画图对比。如果预测曲线和实际偏差大先检查数据是否有异常点再检查N和I0是否合理。参数调整建议如果拟合出的γ对应的感染期明显偏离常识比如超过30天可能是数据里混入了长期阳性病例需要截断或者改用更复杂的模型。β的估计对早期数据敏感可以只用前70%的数据拟合后30%做验证。3.3 模型验证留出法和残差检查拟合完不能只看R0还要做验证。最简单的是留出法用前80%的数据拟合预测后20%计算预测误差。如果误差在可接受范围内比如MAPE小于20%模型可用。另外画残差图如果残差随时间呈现系统性偏差说明模型结构有问题比如需要加入时变β或者干预措施。我一般还会检查S(t)I(t)R(t)是否始终等于N如果数值解出现偏离说明步长太大或者RK4实现有误。这个检查能快速定位低级错误。4. 避坑与排查SIR建模中最容易翻车的五个地方4.1 现象S(t)变成负数曲线出现震荡原因欧拉法步长过大或者RK4实现时dt没有同步更新。SIR方程在I接近0时刚性较强显式方法容易不稳定。解决减小dt到0.01以下或者改用自适应步长的scipy.integrate.solve_ivp。如果必须用固定步长确保RK4的四个斜率都用同一个dt计算。4.2 现象拟合出的R0高得离谱比如大于10原因数据里的现存感染者被低估或者N设置得太小。比如实际人口100万代码里N写成1万β就会被放大。解决核对N的单位和数据单位是否一致。如果数据是每万人中的病例数N就取1。另外检查I0是否用了第一天的数据如果第一天数据是0I0要手动设一个大于0的值。4.3 现象预测曲线比实际晚很多天到达峰值原因γ估计偏小感染期被拉长。常见于数据里包含大量轻症或无症状感染者他们实际移除更快但被计入了长期感染者。解决用累计确诊和累计移除分别拟合或者改用SEIR把潜伏期单独建模。如果坚持用SIR可以对γ加一个下界约束比如gamma 0.05。4.4 现象不同时间段拟合出的β差异巨大原因真实疫情中干预措施戴口罩、限流会改变接触率β不是常数。SIR的常数β假设在长周期内不成立。解决把时间分段每段单独拟合β或者用时变β(t)的扩展模型。如果只是做短期预测用最近两周的数据拟合即可。4.5 现象代码跑得通但结果和文档里的对不上原因文档里的方程可能用了不同的符号约定比如β定义为单位时间接触数而不是传染概率或者R(t)表示累计确诊而不是移除者。解决先用手算一个简单例子验证。取N1000I01β0.3γ0.1手动算前两步再和代码输出对比。这个步骤能省掉后面几小时的排查时间。5. 进阶技巧用SIR做情景分析和干预模拟5.1 把β拆成接触率和传染概率基础SIR里β是一个黑箱参数但实际做干预评估时需要知道改变哪个环节有效。常见做法是把β拆成c人均接触数和p每次接触的传染概率β c * p。这样戴口罩降低p限流降低c可以分别模拟。代码上只需要在sir_deriv里把beta替换成c*p然后对c和p分别扫描。def sir_deriv_intervention(t, y, c, p, gamma, N): S, I, R y beta c * p dS -beta * S * I / N dI beta * S * I / N - gamma * I dR gamma * I return np.array([dS, dI, dR])参数说明c的单位是“接触数/天”p是无量纲概率。如果原来β0.3可以设c10p0.03。干预模拟时把c降到5看峰值如何变化。5.2 用蒙特卡洛评估参数不确定性拟合出的β和γ有置信区间直接拿点估计做预测会低估风险。我一般用蒙特卡洛从拟合的协方差矩阵里采样1000组参数每组跑一次SIR取I(t)的5%和95%分位数作为预测区间。这样报告里能给出“峰值在X到Y之间”的结论比单条曲线可信得多。import numpy as np # 假设res是拟合结果res.hess_inv是逆Hessian近似 params_samples np.random.multivariate_normal(res.x, res.hess_inv, size1000) peak_I [] for beta_s, gamma_s in params_samples: if beta_s 0 or gamma_s 0: continue sol sir_rk4(beta_s, gamma_s, N, I0, R0_init, days, dt0.1) peak_I.append(sol[:, 2].max()) peak_I np.array(peak_I) print(f峰值感染者中位数{np.median(peak_I):.0f}) print(f90%置信区间[{np.percentile(peak_I, 5):.0f}, {np.percentile(peak_I, 95):.0f}])注意res.hess_inv在Nelder-Mead下可能不可用可以改用L-BFGS-B或者用bootstrap重采样。这个技巧在写报告时特别有用能避免“模型预测峰值是1000”这种过于绝对的表述。5.3 一个我常犯的错误早期我做SIR拟合时总想把所有数据点都拟合得很准结果β和γ被少数异常点带偏。后来养成习惯先画散点图肉眼剔除明显离群的天数再拟合。另外SIR模型对早期数据敏感如果第一天只有1个病例I01拟合出的β会偏大。我一般把前3天的数据去掉从第4天开始拟合用前3天估算I0。这个习惯帮我省了很多次返工。希望帮到你。本文还有配套的精品资源点击获取
