- 基础版matlab代码:
clear;clc;
global N; N = 2489.43e4; %地区总人数
T = 1:100; %设置时间观察到第150天
[data,~]=xlsread('data.xlsx');%读取data中的数据放入data中
Ir = data(:,2); %前63天实际感染者数量
z = zeros(length(T)-length(Ir),1); %未知天数数据赋0
Ir = [Ir;z]; %数组扩展
Ir = Ir';
zz = find(Ir==0);
Ir(zz)=NaN;
% Real
% figure
% %plot(T(1:63),Ir);grid on
% plot(T,Ir,'LineWidth',1.5);grid on
% legend('实际感染者')
% xlabel('天数');ylabel('人数')
% hold on
% Estimate
% 定义初始状态各人数
%N = 24894300; %总人数
E = 590; %E类暴露者Exposed,指接触过感染者,但暂无能力传染给其他人的人,对潜伏期长的传染病适用;
%本例中第一天为590个。
I = 401; %I类感病者Infectious,指染上传染病的人,可以传播给S类易感成员,将其变为 E 类或 I 类成员;
%本例中第一天为401个。
S = N - I; %易感者 (Susceptible),指未得病者,但缺乏免疫能力,与感染者接触后容易受到感染;
%假设某区域的人口数为10000,那么第一天的S=N-I=9999
R = 0; %康复者Recovered,指被隔离或因病愈而具有免疫力的人。如免疫期有限,R 类成员可以重新变为 S 类。
%本例中第一天为0个。
r = 20; %感染患者(I)每天接触的易感者数目,本例为20
B = 0.023; %传染系数;由疾病本身的传播能力,人群的防控能力决定,本例设置为0.03
a = 0.07; %潜伏者的发病概率;一般为潜伏期的倒数,本例为1/14
y = 0.1; %恢复系数;一般为病程的倒数,例如流感的病程5天的话,那么它的γ就是1/5,本例设置为1/10=0.1
% T = 1:150; %设置时间观察到第150天
for i = 1:length(T)-1
S(i+1) = S(i) - r*B*S(i)*I(i)/N(1);
E(i+1) = E(i) + r*B*S(i)*I(i)/N(1)-a*E(i);
I(i+1) = I(i) + a*E(i) - y*I(i);
R(i+1) = R(i) + y*I(i);
end
plot(T,S,T,E,T,I,T,R,'LineWidth',1.5);grid on;
xlabel('天数');ylabel('人数')
legend('易感者','潜伏者(暴露者)','感染者','康复者')
% plot(T,I,T,Ir,'LineWidth',1.5);grid on
% xlabel('天数');ylabel('人数')
% legend('预测感染者','实际感染者')
参考来源:https://blog.csdn.net/weixin_39790168/article/details/111008359?spm=1001.2014.3001.5506
参考文献:基于SEIR的新冠肺炎传播模型及拐点预测分析_范如国
- 优化版(增加动态参数)matlab部分代码:
B = 0.5; %初始传染系数
k = 1; %形态参数
L = 22; %衰减速率
for i = 1:length(T)-1
%-----------------------修改部分-------------------------
B = exp(-(T(i)/L)^k*log(2));
S(i+1) = S(i) + B*(S(i)+a*E(i))*I(i)/N(1);
E(i+1) = E(i) + B*(S(i)+a*E(i))*I(i)/N(1)-a*E(i);
%--------------------------------------------------------
%S(i+1) = S(i) - r*B*S(i)*I(i)/N(1);
%E(i+1) = E(i) + r*B*S(i)*I(i)/N(1)-a*E(i);
I(i+1) = I(i) + a*E(i) - y*I(i);
R(i+1) = R(i) + y*I(i);
end
参考文献:基于改进SEIR算法的疫情传播趋势检测方法_燕杨