当前位置：首页 > news >正文

手搓传染病模型（SEIR-拓展）

news 2025/4/29 5:47:12

我们在前文手搓了最基本的SEIR模型

今天我们将根据真实数据拟合出最优参数。

以传染率（β）为例。

话不多说，我们直接开始手搓

% 模型参数
N = 10000; % 总人数
I0 = 10; % 初始感染人数
R0 = 0; % 初始康复人数
E0 = 2 * I0; % 初始暴露人数
S0 = N - I0 - R0 - E0; % 初始易感人数
w = 0.2; % 暴露转化速率
gamma = 1 / 12; % 康复速率
num_days = 100; % 模拟天数
load observed_I.mat % 导入观察病例% 定义拟合选项和初始参数
beta_initial = 0.01;
lb = 0.01;
ub = 1.00;% 使用lsqcurvefit进行非线性最小二乘拟合
options = optimoptions('lsqcurvefit', 'Display', 'iter');
[beta_optimal, resnorm] = lsqcurvefit(@(beta, t) forecast(beta, num_days, N, w, gamma, S0, E0, I0, R0), ...beta_initial,...1: num_days,...observed_I,...lb,...ub,...options);%% 使用最佳的beta进行预测
% x(1):感染人群I， 
% x(2):易感人群S， 
% x(3):康复人群R，
% x(4):暴露人群Edxdt = @(t, x) [w * x(4) - gamma * x(1); % dIdt-beta_optimal * x(2) * x(1) / N; % dSdtgamma * x(1); % dRdtbeta_optimal * x(2) * x(1) / N - w * x(4); % dEdt];[t, y] = ode45(dxdt, 1: num_days, [I0, S0, R0, E0]);hold on
plot(t, y(:, 1));
plot(t, y(:, 2));
plot(t, y(:, 3));
plot(t, y(:, 4));
legend('感染人数I', '易感人数S', '康复人群R', '暴露人群E');

上述代码中的forecast函数见下面。

function I_extracted = forecast(beta, num_days, N, w, gamma, S0, E0, I0, R0)% x(1):感染人群I， % x(2):易感人群S， % x(3):康复人群R，% x(4):暴露人群Edxdt = @(t, x) [w * x(4) - gamma * x(1); % dIdt-beta * x(2) * x(1) / N; % dSdtgamma * x(1); % dRdtbeta * x(2) * x(1) / N - w * x(4); % dEdt];[~, y] = ode45(dxdt, 1: num_days, [I0, S0, R0, E0]);I_extracted = y(:, 1);
end

看看是否可以拟合