Matlab传染病模型实战：SIR/SEIR参数反推与情景推演

简介这份资源面向生物数学、传染病动力学方向的学习者与科研人员围绕Matlab在传染病建模中的应用展开帮助读者理解SIR、SEIR等模型的构建思路并借助代码完成无病平衡点与地方病平衡点的稳定性分析。压缩包共13个文件包含8个m脚本、4个fig图形文件和1份pdf文档整体约2.48MBm文件用于定义模型参数、方程与求解流程fig文件保存了不同参数下的动态仿真结果pdf则提供理论参考。资源中涉及后向分支、E0至E2各平衡点的局部与全局稳定性判定以及基于dde23求解器的延迟微分方程模拟可帮助读者对照代码复现疾病传播规律、验证基本再生数R0对长期行为的影响。目前已有3887人学习下载适合希望将数值计算与传染病防控分析结合的中高级读者参考。1. 用 Matlab 搭传染病模型从 SIR 曲线到参数反推一次讲透很多人第一次在 Matlab 里跑传染病模型都是被一条“S 型曲线”骗进来的——以为把三个微分方程丢进ode45就完事了结果参数一改曲线要么直接躺平要么炸成一条直线完全对不上现实里的感染波峰。传染病模型在 Matlab 里真正难的不是写方程而是把“感染率、恢复率、初始条件”这几个旋钮拧到能复现真实传播节奏的位置并且知道每个旋钮拧过头会发生什么。这篇笔记面向两类人一类是刚接触 SIR/SEIR、想用 Matlab 把模型跑起来的学生和初级工程师另一类是做数据分析、需要快速验证一套传播假设是否站得住脚的从业者。下面从方程怎么落地、参数怎么标定、坑在哪一路写到怎么用模型做情景推演全部给可复现的代码和参数说明。2. SIR 与 SEIR 的方程拆解先搞清楚每个变量代表谁2.1 三个仓室和四个仓室的区别到底在哪SIR 把人群分成三类易感者 S、感染者 I、移除者 R。它的核心假设是“感染后经过一个平均周期就永久移除”适合病程短、康复后免疫明显的场景。SEIR 多了一个暴露者 E也就是已经感染但还没传染性的阶段适合有潜伏期的传播过程。选哪个不是看哪个高级而是看你的数据里有没有“潜伏期”这个可观测的滞后。如果感染曲线相对移除曲线有明显的时间延迟SIR 会把这个延迟硬塞进恢复率里导致参数失真这时候就该上 SEIR。在 Matlab 里SIR 是一个三维状态向量SEIR 是四维。维度一变后面所有参数标定、初值设置、画图逻辑都要跟着改。我一般建议先把 SIR 跑通确认自己对ode45的返回结构、时间向量、参数传递方式都熟了再扩到 SEIR。直接上 SEIR 很容易在“E 的初值给多少”这个问题上卡住因为 E 通常不可直接观测。2.2 把微分方程写成 Matlab 能吃的函数句柄SIR 的标准形式是dS/dt -βSI/NdI/dt βSI/N - γIdR/dt γI。其中 β 是有效接触率γ 是恢复率N 是总人口。下面是一个可以直接保存为sir_model.m的函数function dydt sir_model(t, y, beta, gamma, N) % y(1)S, y(2)I, y(3)R S y(1); I y(2); R y(3); dS -beta * S * I / N; dI beta * S * I / N - gamma * I; dR gamma * I; dydt [dS; dI; dR]; end逻辑说明这里把beta、gamma、N作为额外参数传入而不是写死在函数里是为了后面做参数扫描时不用反复改文件。N放在分母上是把“接触人数”归一化成比例避免人口规模不同导致 β 不可比。注意y是列向量返回的dydt也必须是列向量否则ode45会报维度不一致。参数说明beta的单位是“每人每天有效接触次数”典型范围在 0.1 到 1.0 之间gamma是恢复率的倒数对应平均感染期比如平均 7 天恢复γ≈0.143。N如果设成 100 万初值 I 给 1那早期增长会非常慢因为 S/N 接近 1 但 I 太小这是正常的指数增长起步阶段。2.3 用 ode45 求解并画出三条曲线N 1e6; beta 0.35; gamma 0.14; I0 100; S0 N - I0; R0 0; y0 [S0; I0; R0]; tspan [0 200]; [t, y] ode45((t,y) sir_model(t, y, beta, gamma, N), tspan, y0); S y(:,1); I y(:,2); R y(:,3); plot(t, S, b, t, I, r, t, R, g, LineWidth, 1.5); legend(S,I,R); xlabel(天数); ylabel(人数); title(SIR 模型基本再生数 R0 \beta/\gamma); grid on;逻辑说明ode45是变步长求解器适合这种非刚性方程。匿名函数(t,y) sir_model(...)把额外参数包进去是 Matlab 里传参的标准做法。画图时直接对y的列做索引S 是第一列I 是第二列R 是第三列顺序不能乱。参数说明tspan给 200 天是为了看到完整的感染波峰和回落。如果beta/gamma小于 1感染曲线不会出现明显波峰I 会单调下降这时候要么调大 β要么检查初值 I0 是不是太小导致还没涨起来就结束了。LineWidth设 1.5 是为了在论文或报告里截图清晰不是必须的。3. 参数标定与情景推演让模型对得上数据3.1 基本再生数 R0 怎么算、怎么用R0 β/γ这是判断传播是否扩散的核心指标。R0 1 时感染会增长R0 1 时感染会自然消退。在 Matlab 里可以直接算R0 beta / gamma; fprintf(基本再生数 R0 %.2f\n, R0);逻辑说明这个值不是模型跑出来的而是参数组合出来的。它的意义在于给你一个快速判断如果 R0 只有 1.2那感染波峰不会太高但拖得久如果 R0 到 3 以上波峰又高又陡对医疗资源冲击大。参数说明β 和 γ 的估计通常来自数据拟合。如果手头有每日新增病例可以用lsqcurvefit或fminsearch去最小化模型输出和实际数据的残差。我一般先用fminsearch粗调再用lsqcurvefit精修因为前者不需要工具箱后者对初值敏感但收敛快。3.2 用 fminsearch 做参数反推的最小示例假设你有一组观测到的 I 数据存在I_obs里时间点存在t_obs里function err fit_err(params, t_obs, I_obs, N) beta params(1); gamma params(2); I0 I_obs(1); S0 N - I0; y0 [S0; I0; 0]; [~, y] ode45((t,y) sir_model(t, y, beta, gamma, N), t_obs, y0); I_model y(:,2); err sum((I_model - I_obs).^2); end % 调用 params0 [0.3, 0.1]; params_hat fminsearch((p) fit_err(p, t_obs, I_obs, N), params0);逻辑说明fit_err返回的是残差平方和fminsearch会不断调整params让这个值变小。注意每次迭代都要重新跑ode45所以如果数据点很多、时间跨度很长这个过程会比较慢。可以把t_obs限制在感染上升期和波峰附近减少计算量。参数说明params0是初值猜测给得太离谱会导致收敛到局部极小。我一般会先用网格搜索扫一遍 β 和 γ 的粗略范围比如 β 从 0.1 到 0.6γ 从 0.05 到 0.3找到残差最小的格子再把这个格子中心作为fminsearch的初值。N在整个拟合过程中保持不变不要把它也当成待拟合参数否则会出现参数不可辨识的问题。3.3 情景推演改一个参数看波峰怎么变beta_list [0.2, 0.35, 0.5]; gamma 0.14; figure; hold on; for i 1:length(beta_list) beta beta_list(i); [t, y] ode45((t,y) sir_model(t, y, beta, gamma, N), [0 200], y0); plot(t, y(:,2), LineWidth, 1.5, DisplayName, sprintf(\\beta%.2f, beta)); end legend(show); xlabel(天数); ylabel(感染者人数); title(不同 \beta 下的感染曲线对比); grid on;逻辑说明循环里每次改 β 重新求解把 I 曲线叠在一起。DisplayName配合legend(show)可以自动生成图例不用手动写字符串数组。参数说明β 从 0.2 到 0.5R0 从 1.43 到 3.57波峰高度和出现时间都会明显变化。如果 β 低于 γ也就是 R0 1曲线会直接下降这时候图上几乎看不到波峰不要以为是代码写错了。4. 避坑与排查这五个地方最容易翻车4.1 初值 I0 给太小导致曲线“没反应”现象跑出来的 I 曲线几乎贴着零轴看起来像模型没生效。原因I0 相对于 N 太小比如 N1e6 时 I01早期增长是指数级的但在线性坐标下前几十天都看不出变化。解决要么把 I0 调大到 100 或 1000要么把 y 轴改成对数坐标semilogy这样早期增长就能看清。4.2 把 beta 和 gamma 的单位搞混现象拟合出来的 γ 大于 1或者 β 小得离谱。原因β 的单位是“每人每天”γ 的单位是“每天”如果数据的时间单位是周而不是天γ 会差 7 倍。解决统一时间单位所有参数都按天算。如果原始数据是周报先把时间轴乘以 7 转成天再拟合。4.3 ode45 报“维度不一致”或“输入参数不足”现象运行时报错提示Not enough input arguments或Dimensions of arrays being concatenated are not consistent。原因sir_model的签名和匿名函数里的参数数量对不上或者y0给成了行向量而dydt返回列向量。解决检查sir_model定义里的参数个数确保匿名函数(t,y) sir_model(t, y, beta, gamma, N)传了 5 个参数y0用分号分隔成列向量。4.4 拟合时参数跑到负值现象fminsearch返回的 β 或 γ 是负数物理上无意义。原因fminsearch是无约束优化参数可以取任意实数。解决在fit_err里加惩罚项比如if beta 0 || gamma 0, err 1e10; return; end把负参数直接判为无效。或者改用lsqcurvefit并设置上下界lb和ub。4.5 把 R0 当成固定常数到处用现象不同情景下都用同一个 R0 去解释结论互相矛盾。原因R0 是 β 和 γ 的组合β 会随接触行为变化γ 会随医疗条件变化R0 本身不是常数。解决每次情景推演都重新算 R0并且在报告里写清楚“在当前 β 和 γ 假设下 R0 等于多少”而不是笼统说“R0 是 2.5”。5. 进阶技巧用 SEIR 处理潜伏期用事件函数处理干预5.1 SEIR 的方程和初值设置SEIR 在 SIR 基础上加了 E 仓室dS/dt -βSI/NdE/dt βSI/N - σEdI/dt σE - γIdR/dt γI。σ 是潜伏期转为感染期的速率1/σ 是平均潜伏期。代码结构function dydt seir_model(t, y, beta, sigma, gamma, N) S y(1); E y(2); I y(3); R y(4); dS -beta * S * I / N; dE beta * S * I / N - sigma * E; dI sigma * E - gamma * I; dR gamma * I; dydt [dS; dE; dI; dR]; end逻辑说明E 的初值通常给 0 或者给一个很小的数因为暴露者往往不可观测。如果 E00 且 I00早期 E 会从零开始上升感染曲线会比 SIR 多一个滞后。参数说明σ 的典型值对应潜伏期 3 到 7 天即 σ 在 0.14 到 0.33 之间。β 和 γ 的含义与 SIR 一致。注意 SEIR 的 R0 仍然是 β/γσ 不影响 R0但影响波峰出现的时间。5.2 用事件函数模拟“封控”或“干预”的突然生效Matlab 的ode45支持事件函数可以在某个条件触发时改变参数或停止求解。比如模拟第 30 天开始干预β 从 0.35 降到 0.15function [value, isterminal, direction] intervention_event(t, y) value t - 30; % 第30天触发 isterminal 1; % 停止求解 direction 1; % 只在上穿时触发 end % 第一段0到30天 [t1, y1] ode45((t,y) sir_model(t, y, 0.35, 0.14, N), [0 30], y0, ... odeset(Events, intervention_event)); % 第二段30天之后用第一段的末状态作为初值 y0_new y1(end, :); [t2, y2] ode45((t,y) sir_model(t, y, 0.15, 0.14, N), [30 200], y0_new); % 合并 t_all [t1; t2(2:end)]; y_all [y1; y2(2:end,:)];逻辑说明事件函数返回的value在 t30 时过零isterminal1让求解器停下。然后手动把第一段的末状态取出来作为第二段的初值换一个 β 继续跑。合并时去掉t2的第一个点避免时间重复。参数说明direction1表示只在value从负变正时触发如果写成 0 则双向触发。odeset(Events, ...)是固定写法。第二段的 β 从 0.35 降到 0.15R0 从 2.5 降到 1.07感染曲线会在干预后继续上升一小段然后回落这个“惯性”是 SEIR 和 SIR 都会有的不是 bug。5.3 用热图看 β 和 γ 同时变化时的波峰高度beta_grid 0.1:0.05:0.6; gamma_grid 0.05:0.02:0.3; peak_I zeros(length(beta_grid), length(gamma_grid)); for i 1:length(beta_grid) for j 1:length(gamma_grid) [~, y] ode45((t,y) sir_model(t, y, beta_grid(i), gamma_grid(j), N), [0 300], y0); peak_I(i,j) max(y(:,2)); end end imagesc(gamma_grid, beta_grid, peak_I); colorbar; xlabel(\gamma); ylabel(\beta); title(不同 \beta 和 \gamma 下的感染波峰高度);逻辑说明双重循环遍历参数网格每次求解后取 I 的最大值存进矩阵。imagesc用颜色表示波峰高度横轴 γ纵轴 β。颜色越亮表示波峰越高。参数说明网格步长决定分辨率0.05 和 0.02 是折中再细会明显变慢。tspan给 300 天是为了确保波峰已经出现如果 R0 接近 1波峰可能很晚才来需要再拉长时间。这个热图能帮你快速找到“波峰高度可接受”的参数区间而不是一个个试。5.4 我自己的习惯先固定 N 和 I0再调 β 和 γ跑了这么多轮我最大的教训是不要同时改所有参数。N 和 I0 一旦确定就不要再动否则每次改 β 或 γ 都要重新理解曲线为什么变。我一般先把 N 设成实际人口规模I0 设成 100 左右然后固定 γ 在 0.1 到 0.2 之间只调 β 看波峰怎么变。等 β 的范围摸清了再回头微调 γ。这样每一步只有一个变量在动出了问题能立刻定位。希望帮到你。本文还有配套的精品资源点击获取