用Matlab搭建SIRS传染病模型：从微分方程到参数标定与随机模拟

简介SIRS传染病学模型Matlab完整仿真资源面向从事数学模型仿真、传染病动力学研究的学生与科研人员用于描述易感(Susceptible)、感染(Infectious)、康复(Recovered)、免疫(Immune)四种人群状态的动态转化过程。模型遵循易感个体接触感染者后转变为感染、感染者康复后转为免疫或再次易感的循环机制以四个微分方程分别刻画各状态变化速率综合考虑传染率、康复率、免疫率对疫情流行曲线的影响可作为课程设计、论文复现或防控策略评估的参考实现。压缩包内共4个文件2个m源码即主程序run_SIRS.m与模型函数Fsirs.m结构清晰可直接运行并支持参数调整另含SIRS_model.jpg与1.png示意图便于对照模型结构、代码流程与仿真输出。整体包体仅80KB轻量易用。通过修改传染率、康复率等关键参数可观察感染者峰值、疫情持续时间及群体免疫形成过程帮助深入理解传染病传播机制并为扩展SEIR等更复杂模型打下基础。已有199人浏览学习适合对传染病建模及Matlab实现感兴趣的入门与进阶读者。1. 为什么要自己搭一个SIRS模型而不是直接套用SIR流感季的哨点医院数据曲线总是冲高后回落再过一个季度又抬头。SIR模型里的 R 永远回不到 S解释不了第二波感染。SIRS 把“康复后免疫消退、重新变回易感”这条路径补上用 ρ 参数控制免疫持续时间才能刻画流感的季节性反弹、流感亚型的轮换和肺结核再激活。这篇文章用 Matlab 写一个可读、可改、可拟合哨点数据的 SIRS 模型源码从微分方程讲到数据清洗、参数标定与随机模拟。适合用 Matlab 做传染病动力学、课程设计或疾控数据分析的工程师新手也能按步骤复现。2. SIRS模型的状态方程与Matlab数值求解2.1 三状态转移是SIRS的核心免疫消退由ρ决定SIRS 把人群分成易感 S、感染 I、康复 R 三个仓室假设总人口 N S I R 恒定忽略出生与自然死亡只保留三条事件流易感者以 βSI/N 的速率被感染感染者以 γI 的速率康复康复者以 ρR 的速率失去免疫重新回到易感池。β 是有效接触率γ 是康复率1/γ 是平均感染期ρ 是免疫丧失率1/ρ 是平均免疫持续时间。SIR 模型其实是 SIRS 在 ρ 0 时的特例一旦 ρ 大于 0方程组就没有解析解只能做数值积分。这个方程组反映了传染病动力学的核心机制感染峰的高度由 β 主导峰的宽度由 γ 主导而两波之间的间隔由 ρ 主导。实际操作中很多人直接把 SIR 代码拿来跑 SIRS只加了 ρR 这一项结果发现长期曲线依然单调收敛。原因是 ρ 太小免疫消退在模拟窗口内几乎不起作用。调试时先把 1/ρ 设为观察窗口的一半比如模拟 365 天就设 ρ 1/180立刻能看到第二波的形状。调参顺序应该是先固定 γ再扫 β 确定第一波峰值最后用 ρ 对齐第二波出现的时间点三步分开做不要一次同时动三个参数。2.2 用ode45搭建第一个SIRS求解脚本Matlab 里最稳妥的求解器是 ode45它实现的是变步长四阶 Runge-Kutta 法对 SIRS 这类非刚性问题精度足够速度比 ode15s 快。下面这个最小脚本可以直接保存为 sirs_basic.m 运行。function sirs_basic() % SIRS模型最小求解ode45 匿名函数 beta 0.35; % 有效接触率单位 1/天 gamma 1/7; % 康复率平均感染期7天 rho 1/180; % 免疫丧失率平均免疫维持180天 N 100000; % 总人口 tspan [0 365]; % 模拟一年 y0 [N-10, 10, 0]; % 初始10个感染者 rhs (t,y) [ -beta*y(1)*y(2)/N rho*y(3); beta*y(1)*y(2)/N - gamma*y(2); gamma*y(2) - rho*y(3) ]; [t,y] ode45(rhs, tspan, y0); plot(t, [y(:,1) y(:,2) y(:,3)]*100/N, LineWidth, 1.5); legend(S,I,R,Location,best); xlabel(天); ylabel(占总人口百分比);rhs 是一个匿名函数三行分别对应 S、I、R 三个仓室的导数顺序必须和微分方程的书写顺序一致。y(1)、y(2)、y(3) 是绝对人数而不是百分比所以 rhs 里用 y(1)*y(2)/N 做双线性接触项最后画图时才统一除以 N 再乘 100。如果把方程写成百分比形式双线性项会变成 β 乘两个百分比再乘 N参数的量纲就变了拟合时容易出问题。运行后把 beta 改成 0.5感染峰值会明显提前和增高把 rho 改成 1/60第二波会在当年出现。对比这两个结果就能直观理解 SIRS 与 SIR 的行为差异。注意 ode45 在事件触发时刻不会自动记录状态如果要提取峰值时间需要自己用逻辑判断或 findpeaks。2.3 SIRS的R0推导与长期行为判断对无病平衡点做线性化分析可以得到基本再生数 R0 β/γ这个表达式里没有 ρ因为 ρ 只影响地方病平衡点的位置而不决定疫情是否能够爆发。当 R0 小于等于 1 时感染人数单调下降最后趋近于零当 R0 大于 1 时第一波必然出现之后是否出现第二波取决于 ρ 的大小。R0 越接近 1临界免疫持续期越长此时需要把模拟窗口拉长到 5 年才能看到完整的周期性。% 基于ode45输出判断出现几次感染峰 [~, pk] findpeaks(y(:,2)); fprintf(感染峰数量: %d\n, numel(pk));findpeaks 需要 Signal Processing Toolbox没有工具箱时可以改用循环判断 y(i,2) y(i-1,2) 且 y(i,2) y(i1,2)。判断第二波是否会出现有一个经验式当 1/ρ 小于第一波峰值出现时间到模拟结束时间的间距且 R0 仍大于 1 时模型必然给出第二波。写论文报告时这条结论比直接贴一张双波曲线更有说服力因为它给出了“免疫维持多久会引发下一波”的定量边界。3. SIRS-Matlab源码的模块化组织与哨点数据预处理3.1 源码目录结构与数据文件的组织方式单个脚本跑通以后工程化是下一步。把模型求解、参数拟合、绘图拆成独立函数主脚本只做编排这样换数据、换参数、换绘图样式都不需要动核心方程。一个符合常规的目录结构长这样。sirs_project/ ├── run_all.m ├── src/ │ ├── sirs_rhs.m │ ├── simulate_sirs.m │ ├── fit_sirs_params.m │ └── plot_sirs.m └── data/ ├── raw/ili_raw.csv ├── cleaned/ili_clean.mat └── params/param_sets.csv文件职责输入输出run_all.m主脚本按顺序调用下面三个函数无运行日志sirs_rhs.m微分方程右端函数t, y, 参数结构体dydtsimulate_sirs.m封装 ode45返回时间序列参数结构体, 初始状态t, yfit_sirs_params.m最小二乘标定参数清洗后的观测数据theta_hatplot_sirs.m所有绘图逻辑t, y, 参数图run_all.m 里用 addpath 把 src 和 data 加进搜索路径再依次调用模拟与拟合函数。原始 CSV 和清理后的 MAT 文件分开存放避免每次清洗都覆盖源数据。参数文件用 CSV 保存而不是写在代码里方便记录每一组实验的参数来源和拟合结果。这一套组织方式看似简单实际项目里最常见的故障就是数据不一致某次清洗改了列名另一个脚本还按旧列名读取导致拟合结果全部对不上。统一在 run_all.m 开头用 assert 检查列名和数据长度能省掉半天排查时间。3.2 从ILI%序列到模型输入的3步清洗哨点监测数据通常给的是每周 ILI%也就是流感样病例占门诊就诊量的比例不是 SIRS 需要的感染人数。口径转换加数据清洗一般分三步统一时间戳、处理缺失值和坏点、平滑去噪。function clean_data preprocess_ili(raw_table) % 第1步统一时间戳为datetime排序并去除重复项 t datetime(raw_table.Date, InputFormat, yyyy/MM/dd); [~, idx] sort(t); t t(idx); x raw_table.ILI_percent(idx); % 第2步缺失值前向填充超出合理范围的数据置为NaN x fillmissing(x, previous); x(x 0 | x 20) NaN; % 第3步7点滑动平均消除门诊周末效应 x_smooth movmean(x, 7, omitnan); valid ~isnan(x_smooth); clean_data timetable(t(valid), x_smooth(valid), ... VariableNames, {ILI_percent});fillmissing 用 previous 做前向填充适合短缺口连续缺失超过三个点就应该直接丢弃对应段否则会把跨季节的假方波喂给拟合器。超过 20% 的 ILI% 几乎都是录入错误直接置 NaN。movmean 的 omitnan 选项很关键遇到 NaN 时不会把整个窗口拖成缺失。数据是周度而模型是日度处理方式有两种一是 spline 插值到天二是把模型输出按 7 天聚合后与周报比较。推荐第二种插值会人为制造不存在的日内波动聚合只会丢掉高频细节。3.3 参数网格批量模拟与峰值结果可视化单次模拟只能看一条曲线做情景分析时要扫参数网格。用嵌套循环加上结构体数组保存结果是 Matlab 里最直接的做法。function results run_scenarios(beta_list, rho_list, params) num_b numel(beta_list); num_r numel(rho_list); results struct(beta, cell(num_b,num_r), ... rho, cell(num_b,num_r), ... peak, cell(num_b,num_r), ... t_peak, cell(num_b,num_r)); for i 1:num_b for j 1:num_r p params; p.beta beta_list(i); p.rho rho_list(j); [t, y] simulate_sirs(p); [results(i,j).peak, k] max(y(:,2)); results(i,j).t_peak t(k); results(i,j).beta beta_list(i); results(i,j).rho rho_list(j); end end注意 p params 是结构体的复制必须放在内层循环里重新赋值否则上一次迭代写入的 beta 会泄漏到下一组。这里记录 peak 和 t_peak 两个指标分别代表感染峰值和达峰时间可以区分“把峰压低”和“把峰推迟”这两种干预效果。绘制热力图时用 imagesc 或 heatmap横轴是 beta纵轴是 rho颜色是 peak。实际观察结论通常是 beta 方向颜色变化明显快于 rho这意味着在短周期模拟里接触率比免疫持续期更敏感。数据量较大时把嵌套循环改成 parfor可以把 9 组参数扩展到 9 乘 9 的完整网格。4. SIRS模型参数标定、R0估计与PRCC敏感性分析4.1 用lsqcurvefit把SIRS拟合到观测序列参数标定的目标是最小化模型输出与观测序列之间的均方误差。常规做法是用 Optimization Toolbox 的 lsqcurvefit设置合理的上下界避免参数跑到没有流行病学意义的区间。function fit_sirs_params(data_t, data_y) model (theta, ~) aggregate_output(theta, data_t); theta0 [0.35, 1/7, 1/180]; lb [0.05, 1/14, 1/365]; ub [1.0, 1/2, 1/30]; opts optimoptions(lsqcurvefit, ... Display, iter, ... MaxFunctionEvaluations, 5000, ... FunctionTolerance, 1e-6); theta_hat lsqcurvefit(model, theta0, data_t, data_y, lb, ub, opts); function y_agg aggregate_output(theta, t) beta theta(1); gamma theta(2); rho theta(3); [tm, ym] ode45((t,y) sirs_rhs(t,y,beta,gamma,rho), ... [min(t) max(t)], [99990 10 0]); % 按7天聚合模型输出与周报数据对齐 y_agg movsum(ym(:,2), 7) / 7; y_agg interp1(tm, y_agg, t);theta0 取经验中值保证曲线从第一轮迭代就开始产生可见的波动。上下界要结合疾病常识gamma 的下界 1/14 对应 14 天最长感染期rho 的上界 1/30 对应最短 30 天免疫维持期。aggregate_output 里先做 7 天滑动平均再插值到观测时间点模型输出和数据的口径一致拟合才可能收敛。lsqcurvefit 对初值敏感换一组 theta0 得到不同结果时用多起点拟合随机生成 20 个初值各跑一次取目标函数最小的那组。最后从 theta_hat 直接算 R0 beta_hat / gamma_hat并报告标准误对应的置信区间。4.2 用拉丁超立方加PRCC找关键参数拟合得到的是点估计还需要回答一个问题哪个参数稍微偏一点结果就会大幅变化。偏秩相关系数 PRCC 是传染病模型中常用的全局敏感性指标它把参数和输出分别做秩变换再回归掉其他参数的影响。采样用拉丁超立方保证低样本量下覆盖整个参数空间。n 500; p lhsdesign(n, 3); % 将[0,1]区间映射到参数的物理范围 beta_s p(:,1) * (1.0 - 0.05) 0.05; gamma_s p(:,2) * (1/2 - 1/14) 1/14; rho_s p(:,3) * (1/30 - 1/365) 1/365; peak_out zeros(n, 1); for i 1:n params struct(beta, beta_s(i), ... gamma, gamma_s(i), ... rho, rho_s(i)); [~, y] simulate_sirs(params); peak_out(i) max(y(:,2)); end得到 500 组输入和对应的峰值后对每个参数和输出分别做秩变换然后用线性回归把另外两个参数的影响剔掉残差的相关系数就是 PRCC。实际经验是 beta 的 PRCC 绝对值通常最大因为它在感染项里直接和 S、I 相乘rho 的 PRCC 在一年窗口里比较小但把模拟窗口拉长到五年后明显增大。报告结论时一定要注明时间窗口否则同一个模型会得出相互矛盾的敏感性排序。4.3 SIRS的边界哪些场景必须换模型SIRS 假设均匀混合没有年龄结构也没有空间接触异质性。真实的流感传播有明显的学生与成人接触网络差异SIRS 的 beta 是常数抓不住寒暑假接触矩阵的变化。遇到两类场景建议换模型一是评估学校停课、居家办公这类结构性干预用 SEIR 加接触矩阵二是处理抗体衰减的年龄差异把单室 R 拆成多个免疫等级。SIRS 适合做长期宏观趋势、反复爆发频率和免疫维持时间的数量级判断不适合做逐地精确预测。另外 SIRS 是确定性模型小种群场景下感染可能随机灭绝但 ODE 不会给出灭绝概率这一步必须依赖随机模拟。5. SIRS的随机模拟与Matlab数据对齐的实用技巧5.1 用Gillespie算法模拟随机灭绝当总人口只有几千人时随机波动不能忽略。Gillespie 算法的思路很简单每次迭代根据当前状态计算三个事件的速率用指数分布采样决定下一步发生的时间再按速率比例随机选择事件。下面是一个保留感染历史和易感、康复状态的轻量实现。function [t_hist, I_hist] gillespie_sirs(beta, gamma, rho, S0, I0, R0, tmax) N S0 I0 R0; t 0; S S0; I I0; R R0; t_hist 0; I_hist I0; while t tmax I 0 rates [beta*S*I/N, gamma*I, rho*R]; rsum sum(rates); t t - log(rand()) / rsum; r rand() * rsum; if r rates(1) S S - 1; I I 1; % 感染事件 elseif r rates(1) rates(2) I I - 1; R R 1; % 康复事件 else S S 1; R R - 1; % 免疫消退事件 end t_hist(end1, 1) t; I_hist(end1, 1) I; end这段代码保留了 S、I、R 三个状态的实时变化若只关心感染人数可以只记录 I。批量跑一千次模拟时把 t_hist 和 I_hist 预分配为 1e6 长度的数组再截断否则动态扩展会成为性能瓶颈。同样的参数跑 100 次会有一定比例的轨迹感染峰值很低甚至直接归零这个灭绝概率是确定性模型完全看不到的信息。如果需要输出完整状态用于画多条轨迹可以返回结构体而不是两个数组。5.2 数据对齐与ode45精度取舍一个常被忽略的问题是模型输出时间轴与观测数据的时间基准不一致。周报数据一般落在周一模型输出的时间点是浮点天数直接比较会报时间戳不匹配。用 timetable 的 synchronize 函数统一到周一缺失的时间点自动填补可以避免手工查找对齐错位的数据项。拟合阶段还有一个精度取舍ode45 默认的相对容差是 1e-3画出来曲线略有毛糙但趋势正确做参数标定时先保持默认容差快速搜索参数空间锁定最优参数范围后再把 RelTol 调到 1e-6 精算最终曲线。不要一上来就调小容差批量拟合时每一步多出几倍的计算时间整体进度会慢得让人失去耐心。先用毛糙结果换迭代速度再把干净曲线留给报告。本文还有配套的精品资源点击获取