简介：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(i+1,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, 参数结构体 | dydt |
| simulate_sirs.m | 封装 ode45，返回时间序列 | 参数结构体, 初始状态 | t, y |
| fit_sirs_params.m | 最小二乘标定参数 | 清洗后的观测数据 | theta_hat |
| plot_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(end+1, 1) = t; I_hist(end+1, 1) = I; end这段代码保留了 S、I、R 三个状态的实时变化，若只关心感染人数，可以只记录 I。批量跑一千次模拟时，把 t_hist 和 I_hist 预分配为 1e6 长度的数组再截断，否则动态扩展会成为性能瓶颈。同样的参数跑 100 次，会有一定比例的轨迹感染峰值很低甚至直接归零，这个灭绝概率是确定性模型完全看不到的信息。如果需要输出完整状态用于画多条轨迹，可以返回结构体而不是两个数组。
5.2 数据对齐与ode45精度取舍
一个常被忽略的问题是模型输出时间轴与观测数据的时间基准不一致。周报数据一般落在周一，模型输出的时间点是浮点天数，直接比较会报时间戳不匹配。用 timetable 的 synchronize 函数统一到周一，缺失的时间点自动填补，可以避免手工查找对齐错位的数据项。拟合阶段还有一个精度取舍：ode45 默认的相对容差是 1e-3，画出来曲线略有毛糙但趋势正确；做参数标定时先保持默认容差快速搜索参数空间，锁定最优参数范围后再把 RelTol 调到 1e-6 精算最终曲线。不要一上来就调小容差，批量拟合时每一步多出几倍的计算时间，整体进度会慢得让人失去耐心。先用毛糙结果换迭代速度，再把干净曲线留给报告。
本文还有配套的精品资源，点击获取
