一、簡介
二、源代碼
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);
三、運作結果
四、備注