MATLAB实战传染病动力学建模：从SI、SIS、SIR模型到数值模拟与参数分析

简介本资源是一套面向生物信息学、数学建模与公共卫生研究初学者的MATLAB传染病动力学仿真源码集聚焦SI、SIS、SIR三类经典 compartmental 模型的完整实现解决理论公式到数值模拟的落地难题。压缩包共含多个.m主程序文件如si_sim.m、sis_sim.m、sir_sim.m及配套参数配置脚本涵盖微分方程定义、ode45数值求解、多状态变量时序绘图与关键参数β、γ敏感性分析功能总大小仅75KB轻量易读适合教学演示与课程设计复现。已有4784人学习下载代码结构清晰、注释详实每类模型均提供可直接运行的主函数、状态演化曲线图及参数调整说明帮助读者快速掌握传染病建模核心逻辑、MATLAB动态系统编程规范及结果可视化方法为后续开展疫苗策略仿真或真实疫情数据拟合奠定实践基础。1. 项目概述从零到一用MATLAB玩转传染病动力学建模最近在整理硬盘翻出来一堆以前做数学建模和科研时写的MATLAB代码其中关于传染病模型的这部分我觉得特别有分享价值。SI、SIS、SIR这三个模型可以说是流行病学、信息传播乃至社交网络分析领域的“三原色”几乎所有复杂的传播动力学研究都从它们开始。很多同学在初次接触时往往被一堆微分方程和参数搞得晕头转向或者代码跑通了却不知道结果到底意味着什么。我当年也踩过不少坑比如参数设置不合理导致模型爆炸或者对“基本再生数R0”的理解停留在公式层面。所以我决定把这些源码和背后的思考系统地整理出来。这不仅仅是一份可以“复制粘贴”的代码集更是一次对传染病数学建模核心逻辑的深度拆解。我会带你从最基础的模型假设开始一步步推导微分方程然后用MATLAB实现数值求解和可视化最后再聊聊这些简单模型如何与现实世界对接以及在实际建模竞赛或科研中如何灵活运用和扩展。无论你是正在备战数学建模竞赛的学生还是刚开始接触系统动力学的科研新手相信这份结合了代码与经验的“实战手册”都能让你少走弯路真正理解模型背后的“所以然”。2. 模型基石SI、SIS、SIR的核心思想与方程推导在动手写代码之前我们必须把模型的“地基”打牢。这三个模型描述的都是人群在传染病影响下的状态转移核心是几个关键假设和由此导出的微分方程。2.1 模型基本假设与状态定义所有经典的仓室模型Compartmental Model都始于几个共同假设人口封闭性不考虑出生、死亡非疾病所致和迁移总人口数N恒定。这是一个非常重要的简化让我们能专注于疾病本身的传播动力学。均匀混合人群中个体充分混合任何一个易感者S与任何一个感染者I接触的机会均等。这显然是对复杂社交网络的一种理想化但它是所有基础模型的起点。状态划分根据健康状态将总人口划分为几个互斥的“仓室”。这就是模型名称的由来。S (Susceptible)易感者健康但可能被感染的人群。I (Infectious)感染者已患病且具有传染性的人群。这里我们通常假设感染者一旦进入I状态就立即具有传染性并且传染能力恒定。R (Recovered/Removed)康复者或移除者从感染中恢复并获得永久免疫力或死亡不再参与传播过程。有了这些假设我们就可以用数学语言来描述状态之间如何流动了。流动的速率由关键参数控制。2.2 核心参数β与γ的物理意义模型的行为几乎完全由两个参数决定感染率 β (Beta)这不是一个简单的概率。更准确地说β 有效接触率×每次接触的传染概率。它表示一个感染者单位时间内能成功感染多少个易感者。例如β0.5可以粗略理解为平均每个感染者每天能成功感染0.5个易感者。它的量纲是 1/(时间·人数)但在常微分方程中我们通常处理的是比例所以实际使用的是β/N当总人口N被归一化时或β本身当方程描述的是人数而非比例时。这是新手最容易混淆的地方之一在代码中我们到底该用β还是β/N答案取决于你的方程形式。如果方程是dS/dt -β * S * I那么这里的β已经隐含了“人均”的概念如果方程是dS/dt -β * S * I / N那么这里的β就是绝对感染率。在本文的代码中为直观起见我们采用后一种形式即β表示一个感染者与一个易感者接触并导致感染的概率速率再乘以接触频率。移除率 γ (Gamma)表示感染者单位时间内被移出感染状态康复或死亡的比例。它的倒数 1/γ 具有明确的物理意义平均感染期Duration of Infection。例如γ0.2 /天意味着平均感染期是5天。γ越大病人好得或离开得越快。一个极其重要的衍生参数是基本再生数 R0。它的定义是在完全易感的人群中一个感染者在其整个传染期内平均能感染的人数。计算公式很简单R0 β / γ。R0是决定疫情走向的“阈值”R0 1疾病可能流行R0 1疾病会自然消亡。理解R0是理解所有后续模拟结果的关键。2.3 微分方程推导从文字描述到数学公式现在我们把文字描述变成微分方程。设S(t), I(t), R(t)分别表示t时刻易感者、感染者和移除者的人数总人口N S I RSIR模型或 N S ISI/SIS模型。SI模型最简单只有感染没有恢复和免疫。假设感染者终身具有传染性且不恢复如某些疱疹病毒的理想化模型。流程S - I方程dS/dt -β * S * I / NdI/dt β * S * I / N解读易感者减少的速率-dS/dt等于新感染者的产生速率dI/dt这个速率与当前易感者数量S、感染者数量I成正比与总人口N成反比因为接触机会与人口密度有关这里其实是均匀混合假设的数学体现β/N确保了当S和I是人数时转移速率合理。SIS模型感染后可恢复但恢复后不产生免疫力会再次变为易感者如普通感冒、细菌性性病。假设康复后立即再次易感。流程S - I - S方程dS/dt -β * S * I / N γ * IdI/dt β * S * I / N - γ * I解读比SI模型多了一项γ * I表示感染者以速率γ恢复并重新加入易感者行列。这个模型可能达到一个地方性流行平衡点Endemic Equilibrium而不是所有人最终都被感染。SIR模型经典且最重要的模型感染后获得永久免疫力如天花、麻疹、水痘。假设康复后获得永久免疫不再参与传播。流程S - I - R方程dS/dt -β * S * I / NdI/dt β * S * I / N - γ * IdR/dt γ * I解读这是Kermack-McKendrick模型的基石。易感者不断转化为感染者感染者以速率γ转化为移除者。疫情最终会停止但并非所有易感者都会被感染会有一部分人S(∞)从未感染。最终感染规模取决于R0和初始条件。注意在推导和编程时务必检查方程是否满足总人口守恒。对SIR模型d(SIR)/dt 0对SIS模型d(SI)/dt 0。这是一个快速验证方程书写是否正确的好方法。3. MATLAB实战从方程到动态模拟理论清晰后我们进入实战环节。用MATLAB求解这类常微分方程组ODEode45函数是我们的首选利器。它是一种自适应步长的Runge-Kutta方法对于非刚性的传染病模型精度和效率都很不错。3.1 环境准备与代码框架搭建首先我们明确代码结构。一个好的习惯是为每个模型编写独立的函数文件然后在主脚本中调用并绘图对比。这样代码清晰易于管理和调试。% 主脚本 main_epidemic_models.m % 清空环境 clear; close all; clc; % 设置全局参数也可以作为函数参数传递这里用全局变量方便演示 global beta gamma N beta 0.3; % 感染率 gamma 0.1; % 恢复率 N 1000; % 总人口 % 计算基本再生数 R0 R0 beta / gamma; fprintf(基本再生数 R0 %.2f\n, R0); % 初始条件 [S0, I0, R0] I0 1; % 初始感染者 S0 N - I0; % 初始易感者 R0_init 0; % 初始康复者 % 时间跨度单位天 tspan [0, 150]; % 分别调用不同模型的求解函数 [t_SI, Y_SI] solve_SI(tspan, [S0, I0]); [t_SIS, Y_SIS] solve_SIS(tspan, [S0, I0]); [t_SIR, Y_SIR] solve_SIR(tspan, [S0, I0, R0_init]); % 绘制结果 plot_models(t_SI, Y_SI, t_SIS, Y_SIS, t_SIR, Y_SIR);这个主脚本设置了公共参数计算了关键的R0并规划了模拟的时间范围。接下来我们实现三个核心的模型函数。3.2 SI模型实现疫情无可避免的蔓延SI模型描述了最悲观的情况一旦感染永不康复最终所有人都会被感染。function [t, Y] solve_SI(tspan, y0) % 求解SI模型 % 输入 % tspan - 时间范围如 [0, 100] % y0 - 初始条件向量 [S0, I0] % 输出 % t - 时间点向量 % Y - 解矩阵第一列是S第二列是I global beta N % 定义SI模型的微分方程组 function dydt ode_SI(t, y) S y(1); I y(2); dSdt -beta * S * I / N; dIdt beta * S * I / N; dydt [dSdt; dIdt]; end % 使用ode45求解 [t, Y] ode45(ode_SI, tspan, y0); end关键点解析我们使用嵌套函数ode_SI来定义方程这样可以方便地访问全局参数beta和N。方程dSdt -beta * S * I / N中的beta/N是关键。这里beta可以理解为“一个感染者单位时间内有效接触的人数”而(S/N)是接触对象为易感者的概率。因此beta * I * (S/N)就是新感染发生的总速率。ode45会返回一系列时间点t和对应的状态矩阵Y。Y的每一行对应一个时间点每一列对应一个状态变量。3.3 SIS模型实现地方性流行的平衡SIS模型引入了恢复机制但无免疫因此疾病可能持续存在于人群中。function [t, Y] solve_SIS(tspan, y0) % 求解SIS模型 global beta gamma N function dydt ode_SIS(t, y) S y(1); I y(2); dSdt -beta * S * I / N gamma * I; % 减少的易感者 恢复的感染者 dIdt beta * S * I / N - gamma * I; % 新增的感染者 - 恢复的感染者 dydt [dSdt; dIdt]; end [t, Y] ode45(ode_SIS, tspan, y0); end与SI模型的区别在dSdt中多了一项 gamma * I代表康复者重新变为易感者。在dIdt中相应地从新增项里减去了恢复项- gamma * I。这个系统存在一个非零的平衡点。当dIdt 0时可以解得平衡时的感染者比例I* / N 1 - 1/R0当 R0 1。这意味着只要R0大于1无论初始有多少感染者疾病都会稳定在一个固定的流行水平而不是感染所有人。3.4 SIR模型实现经典的疫情浪潮与群体免疫SIR模型是分析疫情爆发、高峰和终结的最有力工具。function [t, Y] solve_SIR(tspan, y0) % 求解SIR模型 global beta gamma N function dydt ode_SIR(t, y) S y(1); I y(2); R y(3); dSdt -beta * S * I / N; dIdt beta * S * I / N - gamma * I; dRdt gamma * I; dydt [dSdt; dIdt; dRdt]; end [t, Y] ode45(ode_SIR, tspan, y0); endSIR模型的核心洞察疫情高峰感染者数量I(t)会先上升后下降。高峰出现在dI/dt 0的时刻此时S N * (γ / β) N / R0。也就是说当易感者人口减少到总人口的1/R0时疫情达到顶峰。最终规模疫情不会感染所有人。当I - 0时易感者比例S(∞)/N是一个大于0的数。最终感染规模R(∞)可以通过求解一个超越方程得到它总是小于总人口。这直观地解释了“群体免疫”的概念当足够多的人通过感染或接种疫苗获得免疫进入R仓室后病毒的传播链就会被打断。相平面分析SIR模型可以简化为dI/dS -1 (N / (R0 * S))积分后得到I S - (N/R0)*ln(S) constant。这个关系式决定了(S, I)相平面上的轨迹是理论分析的有力工具。3.5 可视化与结果分析让数据说话代码的最后一环是将结果清晰呈现。我们编写一个绘图函数来对比三个模型。function plot_models(t_SI, Y_SI, t_SIS, Y_SIS, t_SIR, Y_SIR) global N figure(Position, [100, 100, 1200, 800]); % 设置大图窗 % 绘制SI模型 subplot(2, 3, 1); plot(t_SI, Y_SI(:,1), b-, LineWidth, 2); hold on; % S plot(t_SI, Y_SI(:,2), r-, LineWidth, 2); % I xlabel(时间 (天)); ylabel(人数); title(SI模型 - 最终全部感染); legend(易感者 S, 感染者 I, Location, best); grid on; % 绘制SIS模型 subplot(2, 3, 2); plot(t_SIS, Y_SIS(:,1), b-, LineWidth, 2); hold on; plot(t_SIS, Y_SIS(:,2), r-, LineWidth, 2); xlabel(时间 (天)); ylabel(人数); title(sprintf(SIS模型 - 地方性流行 (平衡点 I*%.0f), N*(1-1/(beta/gamma)))); legend(易感者 S, 感染者 I, Location, best); grid on; % 绘制SIR模型 subplot(2, 3, 3); plot(t_SIR, Y_SIR(:,1), b-, LineWidth, 2); hold on; plot(t_SIR, Y_SIR(:,2), r-, LineWidth, 2); plot(t_SIR, Y_SIR(:,3), g-, LineWidth, 2); xlabel(时间 (天)); ylabel(人数); title(SIR模型 - 疫情爆发与终结); legend(易感者 S, 感染者 I, 康复者 R, Location, best); grid on; % 绘制SIR模型的相平面图 (S-I图) subplot(2, 3, 4); plot(Y_SIR(:,1), Y_SIR(:,2), k-, LineWidth, 1.5); xlabel(易感者 S); ylabel(感染者 I); title(SIR模型相平面 (S-I)); grid on; % 标记初始点和方向 hold on; plot(Y_SIR(1,1), Y_SIR(1,2), ro, MarkerSize, 10, MarkerFaceColor, r); quiver(Y_SIR(1:10:end,1), Y_SIR(1:10:end,2), ... gradient(Y_SIR(1:10:end,1)), gradient(Y_SIR(1:10:end,2)), ... 0.5, Color, [0.5 0.5 0.5]); % 简化方向场示意 % 绘制不同R0下的最终感染规模对比SIR subplot(2, 3, [5,6]); R0_range 0.5:0.1:5; final_R_frac zeros(size(R0_range)); for i 1:length(R0_range) R0_val R0_range(i); % 对于给定的R0求解最终康复者比例近似解通过模拟 beta_temp R0_val * gamma; % 固定gamma改变beta % 这里简化计算使用近似公式最终感染规模比例 ≈ 1 - exp(-R0 * 最终规模比例) % 采用迭代法求解超越方程 s_inf exp(-R0*(1-s_inf))其中s_inf S(∞)/N s_inf 1; % 初始猜测 for iter 1:100 s_new exp(-R0_val * (1 - s_inf)); if abs(s_new - s_inf) 1e-6 break; end s_inf s_new; end final_R_frac(i) 1 - s_inf; % R(∞)/N end plot(R0_range, final_R_frac, m-, LineWidth, 2); xlabel(基本再生数 R_0); ylabel(最终感染人口比例 R(\infty)/N); title(SIR模型最终感染规模 vs R_0); grid on; hold on; plot([1,1], [0,1], r--); % 标记R01的阈值线 text(1.1, 0.1, R_01 (流行阈值), Color, r); plot(R0_range, 1 - 1./R0_range, b--); % 群体免疫阈值线 legend(最终感染比例, 流行阈值, 群体免疫阈值 (1-1/R_0), Location, southeast); sgtitle(传染病基础模型动力学对比 (SI, SIS, SIR)); end运行主脚本你会得到一张综合性的分析图。从图中可以清晰地看到SI模型易感者单调递减至0感染者单调递增至总人口N。SIS模型易感者和感染者数量震荡后趋于一个非零的稳定平衡直观展示了地方性流行。SIR模型经典的疫情曲线。感染者先升后降易感者持续减少康复者持续增加直至疫情结束。相平面图展示了状态空间的演化轨迹。最右边的子图则定量揭示了R0对最终疫情规模的巨大影响当R0略大于1时最终感染比例可能并不高但当R0达到3时最终可能感染超过90%的人口。图中蓝色的虚线1-1/R0就是著名的“群体免疫阈值”当免疫人口比例超过这个值时疫情就会开始衰退。4. 参数敏感性分析与模型扩展初探跑通基础模型只是第一步。在真实的研究或竞赛中我们更需要探究模型的行为如何随参数变化以及如何扩展模型以描述更复杂的现象。4.1 关键参数β和γ的影响β感染率和γ恢复率是模型的“方向盘”和“刹车”。我们可以通过简单的参数扫描来观察它们的影响。% 参数敏感性分析脚本 sensitivity_analysis.m clear; close all; clc; N 1000; I0 1; S0 N - I0; R0_init 0; tspan [0, 200]; % 案例1固定gamma改变beta即改变R0 gamma_fixed 0.1; beta_range [0.05, 0.2, 0.4]; % 对应 R0 0.5, 2, 4 figure; for i 1:length(beta_range) beta beta_range(i); [t, Y] solve_SIR(tspan, [S0, I0, R0_init]); subplot(2,2,1); plot(t, Y(:,2), LineWidth, 1.5, DisplayName, sprintf(\\beta%.2f, R0%.1f, beta, beta/gamma_fixed)); hold on; end subplot(2,2,1); xlabel(时间 (天)); ylabel(感染者 I); title(不同感染率 \beta (固定 \gamma0.1) 对疫情曲线的影响); legend(show); grid on; % 案例2固定R0改变gamma即同时改变beta R0_fixed 2; gamma_range [0.05, 0.1, 0.2]; for i 1:length(gamma_range) gamma gamma_range(i); beta R0_fixed * gamma; [t, Y] solve_SIR(tspan, [S0, I0, R0_init]); subplot(2,2,2); plot(t, Y(:,2), LineWidth, 1.5, DisplayName, sprintf(\\gamma%.2f, \\beta%.2f, gamma, beta)); hold on; end subplot(2,2,2); xlabel(时间 (天)); ylabel(感染者 I); title(不同恢复率 \gamma (固定 R02) 对疫情曲线的影响); legend(show); grid on;分析结果改变β固定γβ越大R0越大疫情峰值越高、来得越早最终感染规模也越大。R01时疾病无法流行曲线几乎贴着0轴。改变γ固定R0为了保持R0不变β需要成比例变化。此时γ越大意味着病程越短疫情峰值越低但爆发和结束得越快曲线更“瘦高”。这解释了为什么缩短感染期如通过有效治疗是控制疫情的重要手段。4.2 模型扩展方向与思路基础SIR模型是骨架现实情况需要添加更多细节。以下是几个常见的扩展方向也是数学建模竞赛中常考的考点SEIR模型在S和I之间增加一个“潜伏期Exposed”仓室E。感染者被感染后不会立即具有传染性而是先进入潜伏期。方程变为S - E - I - R。这需要引入一个新的参数潜伏期转化为传染期的速率σσ的倒数即平均潜伏期。这个模型能更好地描述像COVID-19、流感这样有显著潜伏期的疾病。带出生死亡的SIR模型考虑自然出生率Λ和死亡率μ。方程中会增加ΛN到dS/dt并对所有仓室增加-μ*仓室的项。这可以研究疾病在人口中的长期存在性是消亡还是成为地方病。年龄结构模型将人口按年龄分组不同年龄组的接触率β、恢复率γ甚至死亡率都可能不同。这需要用一个接触矩阵来描述不同组间的接触模式模型会从常微分方程ODE变为偏微分方程PDE或高维ODE系统。这是研究疫苗接种策略优先给哪个年龄组接种的关键。空间异质性模型用元胞自动机Cellular Automata或网络模型Network Model代替均匀混合假设。每个个体是网络中的一个节点疾病沿边传播。这能研究社交距离、社区结构对传播的影响。随机性模型基础ODE是确定性的但现实传播充满随机性。可以用随机微分方程SDE或基于Gillespie算法的随机模拟来研究小规模初始疫情灭绝的概率、疫情规模的波动等。实操心得在数学建模竞赛中不要一开始就追求最复杂的模型。从SIR或SEIR这样的基础模型出发先跑出结果画出图完成基本分析。然后再根据题目要求有选择地引入一两个扩展比如如果题目提到“潜伏期”就加E仓室如果提到“不同人群”就考虑年龄结构或接触矩阵。这样既能保证有扎实的基础分又能体现你的建模深度。代码实现上建议先写好基础模型的通用求解框架扩展时只需修改定义微分方程组的函数即可绘图和分析代码可以复用。5. 常见问题与调试技巧实录在实际编写和运行这些模型时你肯定会遇到各种问题。下面是我总结的一些“坑”和解决方法。5.1 数值求解器报错或不稳定问题使用ode45求解时可能出现“积分容差无法满足”的错误或者结果出现负值、数值爆炸如感染者人数超过总人口。原因与排查参数值不合理这是最常见的原因。β和γ通常应在0到1之间以“每天”为单位。如果你不小心设置了β10意味着一个感染者每天能感染10个人疫情会瞬间爆炸导致求解器步长过小或失败。务必检查参数的数量级。一个经验法则是R0通常在1-10之间γ的倒数平均感染期对于流感是3-7天对于麻疹是7-14天。初始条件不合理确保初始各仓室人数之和等于总人口N。如果S0 I0 R0 N模型会出问题。刚性问题如果参数差异极大例如β很大而γ很小方程可能变成“刚性”的。ode45适用于非刚性方程对于刚性问题会效率低下或失败。可以尝试使用适用于刚性问题的求解器如ode15s或ode23s。% 将 ode45 替换为 ode15s options odeset(RelTol,1e-6, AbsTol,1e-9); % 可以调整容差 [t, Y] ode15s(ode_SIR, tspan, y0, options);方程书写错误这是最隐蔽的错误。务必用总人口守恒来验证对于SIR计算d(SIR)/dt是否恒为0。可以在ode函数末尾加一句assert(abs(sum(dydt)) 1e-10, 方程不守恒)来辅助调试正式运行时注释掉。5.2 结果与理论预期不符问题模拟出的曲线形状奇怪或者最终值不符合理论计算例如SIS模型没有稳定在预测的平衡点。排查步骤检查R0首先打印出你使用的R0值。用公式R0 β / γ手动验算一下。验证平衡点对于SIS模型计算理论平衡点I_star N * (1 - 1/R0)。在你的模拟结果中观察疫情后期I(t)是否在I_star附近小幅波动由于数值误差可能不会完全相等。如果R01I_star为负理论上疾病会消亡模拟中I(t)应趋近于0。延长模拟时间有时疫情发展较慢你设置的tspan太短系统还没有达到稳定状态。尝试将模拟时间延长例如到500天或1000天。检查归一化这是另一个常见错误。如果你的微分方程是针对人口比例sS/N, iI/N, rR/N写的那么方程形式会变为ds/dt -β * s * i其中β的意义与之前不同。确保你代码中的方程形式、参数意义和初始条件是人数还是比例是自洽的。我强烈建议在初学时始终使用人数而非比例来编写方程并在所有涉及人口的地方显式地除以N这样物理意义最清晰不易出错。5.3 性能优化与代码整洁问题当需要做大量参数扫描或运行复杂扩展模型时代码速度慢或者结构混乱难以维护。优化技巧向量化参数扫描避免在循环内频繁调用ode45。可以尝试一次求解多个初始条件或参数但这需要改写ODE函数。更简单实用的方法是使用parfor并行循环如果拥有并行计算工具箱。beta_list linspace(0.1, 0.5, 20); peak_I zeros(size(beta_list)); parfor idx 1:length(beta_list) % 改为 parfor 进行并行计算 beta beta_list(idx); [t, Y] solve_SIR(tspan, [S0, I0, R0_init]); % solve_SIR函数内部需能接收参数 peak_I(idx) max(Y(:,2)); end使用函数句柄传递参数避免使用global全局变量。更优雅的方式是将参数作为额外参数传递给ODE函数。function dydt ode_SIR_param(t, y, beta, gamma, N) S y(1); I y(2); R y(3); dSdt -beta * S * I / N; dIdt beta * S * I / N - gamma * I; dRdt gamma * I; dydt [dSdt; dIdt; dRdt]; end % 调用时 [t, Y] ode45((t,y) ode_SIR_param(t,y,beta,gamma,N), tspan, y0);模块化设计就像本文的示例一样将模型求解、参数设置、绘图分析分离成不同的函数或脚本。一个models文件夹存放solve_SI.m,solve_SIS.m,solve_SIR.m一个utils文件夹存放plot_results.m,calc_R0.m等。主脚本只负责协调。这极大提高了代码的可读性和可复用性。5.4 从模拟到现实参数估计的挑战问题模型很好但β和γ这些参数从哪里来如何让模型拟合真实数据思路这是传染病建模从“玩具”走向“实用”的关键一步。通常有两种途径从文献中获取对于已知疾病其平均感染期1/γ和基本再生数R0的估计值可以在学术论文或权威卫生机构如WHO、CDC的报告中找到。然后利用β R0 * γ推算感染率。从数据中拟合如果你有真实的时间序列数据如每日新增感染人数你可以定义模型输出如每日新增病例β * S * I / N与真实数据之间的误差如最小二乘误差然后使用优化算法如MATLAB的fminsearch,lsqcurvefit来寻找最优的β和γ以及可能的初始条件I0。这是一个反问题通常比较复杂且结果对数据质量和模型假设非常敏感。踩坑记录我曾试图用早期非常粗糙的COVID-19数据拟合SIR模型结果发现拟合出的R0波动极大。原因在于早期检测能力不足数据存在严重的低估和滞后。教训是模型结果的可靠性极度依赖于输入数据的质量。在建模论文中必须对数据来源、局限性和可能存在的偏差进行讨论。对于SIR模型另一个常见问题是它假设康复后终身免疫且传染力恒定这与许多疾病的实际情况如抗体衰减、症状前传染不符。这时就需要选择或构建更合适的模型如SEIR, SIRS等。理解模型的局限性和知道如何用它同样重要。本文还有配套的精品资源点击获取