Prediction Model: Epidemic analysis and prediction based on MATLAB SEIR model

clear; clc; [data,~]=xlsread('data.xls'); % read data into data % declare the global variable global g; % latent incidence global t; % Isolation rate global MU; % cure probability global S; % Number of susceptible people % Number of sleeper global Ii; % Number of cases but not isolated global Iu; % Number of cases and quarantines global D; % death toll global R; % exit now_confirmed=data(45, :); % Confirmed sofar_confirmed=data(46, :); % death=data(47, :); % mobile death cured=data(48, :); % cumulative cure N=5917e4; % Total population [testData,~]=xlsread('test.xlsx'); 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'.'Cure rate');
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('real'.'fit');

%estimate E----------------------------------------------------------------
figure('name'.'Sleeper count');
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('real'.'fit');
E=f(coeff,1:40); E=E(1:end2 -);


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


%--------------------------------------------------------------------------
D=D(1:end2 -);
R=R(1:end2 -);
Iu=Iu(1:end2 -);
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)]; % define initial state [x fval]=fmincon('ff'[0.05 0.05 1e-7 1e-7], [], [], [], [], [0 0 0 0], [5 5 0.1 0.1]); % To find the minimum value mu=polyval(p,1:100);
y=SEIR_discrete(start,t,g,mu,x(1),x(2),x(3),x(4),54);
Copy the code

Version: 2014 a