天天看點

【預測模型】基于matlab SEIR模型疫情分析預測【含Matlab源碼 666期】

一、簡介

二、源代碼

clear;clc;
[data,~]=xlsread('data.xls');%讀取data中的資料放入data中
%聲明全局變量
global g;           %潛伏着發病率
global t;           %發病者隔離率
global mu;          %治愈機率
global S;           %易感者人數
global E;           %潛伏者數量
global Ii;          %發病但未隔離數量
global Iu;          %發病且隔離數量
global D;           %死亡人數
global R;           %退出者

now_confirmed=data(45,:);       %現存确診
sofar_confirmed=data(46,:);     %累計确診
death=data(47,:);               %累計死亡
cured=data(48,:);               %累計治愈
N=5917e4;                       %總人口數
[testdata,~]=xlsread('test.xlsx');%讀取test中的資料放入testdata中
Iu_new=testdata(1,:);
R_new=testdata(4,:);
D_new=testdata(3,:);
%---------------------------------------------------


D=[death(25:end) D_new];
R=[cured(25:end) R_new];
Iu=[now_confirmed(25:end) Iu_new];
sofar_confirmed=[sofar_confirmed(25:end) testdata(2,:)];
figure
plot(R,'g');
hold on;
plot(D,'r');
hold on;
plot(Iu,'b');


%estimate mu---------------------------------------------------------------
figure ('name','治愈率');
plot(diff(R)./Iu(2:end),'r');hold on;
p=polyfit(1:length((diff(R))./Iu(2:end)),diff(R)./Iu(2:end),1);
plot(1:length(diff(R)./Iu(2:end)),polyval(p,1:length(diff(R)./Iu(2:end))),'b');
mu=polyval(p,1:length(diff(R)./Iu(2:end)));
legend('實際','拟合');

%estimate E----------------------------------------------------------------
figure('name','潛伏者數量');
E=[];
for i=1:28
    E(end+1)=sofar_confirmed(i+12)-sofar_confirmed(i);
end
plot(E,'b');hold on
f=@(x,xdata)(x(1)*exp(-1*(x(2)*xdata+x(3))));
coeff=lsqcurvefit(f,[16000,1,0],1:length(E),E);
plot(f(coeff,1:length(E)),'r');
title('E estiamte');
legend('實際','拟合');
E=f(coeff,1:40);E=E(1:end-2);


%estimate Ii---------------------------------------------------------------
Ii=[];
for i=1:38
    Ii(end+1)=sofar_confirmed(i+2)-sofar_confirmed(i);
end


%--------------------------------------------------------------------------
D=D(1:end-2);
R=R(1:end-2);
Iu=Iu(1:end-2);
S=N-D-R-Ii-Iu-E;
t=1;
g=0.1;
global start;
start=[S(1),E(1),Ii(1),Iu(1),R(1),D(1)];            %定義初始狀态
[x fval]=fmincon('ff',[0.05 0.05 1e-7 1e-7],[],[],[],[],[0 0 0 0],[5 5 0.1 0.1]);       %求解非線性多元函數最小值
mu=polyval(p,1:100);
y=SEIR_discrete(start,t,g,mu,x(1),x(2),x(3),x(4),54);
           

三、運作結果

四、備注