武汉加油——传染病模型拟合

武汉加油——传染病模型拟合

最近看了一篇很棒的文章，解决了我许多关于传染病模型困惑。有兴趣请看：
http://i.dataguru.cn/mportal.php?mod=view&aid=15351
因此想对这篇文章的做法复现下，用matlab完成。
思路：使用SIR模型，对参数进行合理拟合，使最终SIR模型仿真尽量接近真实结果，使用数据是武汉地区新冠肺炎的数据。
参数拟合使用回溯传播模型，详细原理见文章，代码将在以下一并展示。
仿真结果：

SIR模型：

gamma=1/C=1/14,gamma是治愈率,C是治愈周期。
beta=k*b,k是每个病人每天接触的人，b是传染概率,beta是每个病人每天有效接触的人
左图是用回溯传播模型拟合参数beta，右图使用公布数据拟合参数beta。k=5。

下图展示的是不同方法拟合参数b，同时k分别为1或1.5的仿真结果,表示政府管控的力度强弱不同的情况。

据文章分析，第一个图的情况比较贴近早期情况。

代码如下：

%%%%%%%%%%%%%%%%参数估计###########33

d=0.03;%死亡率
C=14;
global N gamma
gamma=1/C;
N=13000000;

t0=0;
t1=35;
t2=41;
I0=0;
I1=log(1182);
I2=log(2758);

global k;
k=1;
global k1;
k1=1.5;
t=[t0,t1,t2];
I=[I0,I1,I2];
p=polyfit(t,I,1);
b=(p(1)+gamma)/k;
global beta beta2;
beta=b*k;
beta2=-b*k1;

t=42:46;
t=[0,t];
I=[0,log(198),log(218),log(320),log(478),log(639)];
p1=polyfit(t,I,1);
b1=(p1(1)+gamma)/k;
global beta1 beta3;
beta1=b1*k;%公布
beta3=b1*k1;%公布
%%%%%%%%%SIR模型##############
options = odeset('MaxStep',0.001,'NonNegative',1:3);
%[T,Y] = ode45(@func_SIR, [0,100], [59169802 198 0], options);%[T,Y] = solver(odefun,tspan,y0)。
[T,Y] = ode45(@func_SIR, [0 230], [N-1;1;0], options);

D=zeros(length(Y),1);
for i=C+2:length(Y)
    D(i)=d*(N-Y(i-C-1,1));
end
    
subplot(2,2,1);
plot(T,Y(:,1),'r','LineWidth',2);
hold on;
plot(T,Y(:,2),'k','LineWidth',2);
hold on;
plot(T,Y(:,3),'g','LineWidth',2);
hold on;
plot(T,D,'b','LineWidth',2);

xlabel('t');
ylabel('people');
title('SIR,k=1')
legend({'S','I','R','D'},'Location','SouthEast');%'southeastoutside'

%%%%%%%%%SIR模型公布人少##############
options = odeset('MaxStep',0.001,'NonNegative',1:3);
%[T,Y] = ode45(@func_SIR, [0,100], [59169802 198 0], options);%[T,Y] = solver(odefun,tspan,y0)。
[T,Y] = ode45(@func_SIR1, [0 230], [N-1;1;0], options);

D=zeros(length(Y),1);
for i=C+2:length(Y)
    D(i)=d*(N-Y(i-C-1,1));
end
    
subplot(2,2,2);
plot(T,Y(:,1),'r','LineWidth',2);
hold on;
plot(T,Y(:,2),'k','LineWidth',2);
hold on;
plot(T,Y(:,3),'g','LineWidth',2);
hold on;
plot(T,D,'b','LineWidth',2);

xlabel('t');
ylabel('people');
title('SIR公布k=1')
legend({'S','I','R','D'},'Location','SouthEast');%'southeastoutside'

%%%%%%%%%SIR模型##############
options = odeset('MaxStep',0.001,'NonNegative',1:3);
%[T,Y] = ode45(@func_SIR, [0,100], [59169802 198 0], options);%[T,Y] = solver(odefun,tspan,y0)。
[T,Y] = ode45(@func_SIR2, [0 230], [N-1;1;0], options);

D=zeros(length(Y),1);
for i=C+2:length(Y)
    D(i)=d*(N-Y(i-C-1,1));
end
    
subplot(2,2,3);
plot(T,Y(:,1),'r','LineWidth',2);
hold on;
plot(T,Y(:,2),'k','LineWidth',2);
hold on;
plot(T,Y(:,3),'g','LineWidth',2);
hold on;
plot(T,D,'b','LineWidth',2);

xlabel('t');
ylabel('people');
title('SIR,k=1.5')
legend({'S','I','R','D'},'Location','SouthEast');%'southeastoutside'

%%%%%%%%%SIR模型##############
options = odeset('MaxStep',0.001,'NonNegative',1:3);
%[T,Y] = ode45(@func_SIR, [0,100], [59169802 198 0], options);%[T,Y] = solver(odefun,tspan,y0)。
[T,Y] = ode45(@func_SIR3, [0 230], [N-1;1;0], options);

D=zeros(length(Y),1);
for i=C+2:length(Y)
    D(i)=d*(N-Y(i-C-1,1));
end
    
subplot(2,2,4);
plot(T,Y(:,1),'r','LineWidth',2);
hold on;
plot(T,Y(:,2),'k','LineWidth',2);
hold on;
plot(T,Y(:,3),'g','LineWidth',2);
hold on;
plot(T,D,'b','LineWidth',2);

xlabel('t');
ylabel('people');
title('SIR公布,k=1.5')
legend({'S','I','R','D'},'Location','SouthEast');%'southeastoutside'


function dy=func_SIR(t,y)
global N gamma beta
%beta=77/(198*(3^2-1));%设接触感染者发病概率（1.19,67+10=77）
S = y(1);
I = y(2);
R = y(3);
%x=beta0*(cos(2*pi*tt)*beta1+1);
dS = -beta*I*S/N;
dI =  beta*I*S/N-gamma*I;
dR = gamma*I;
dy=[dS;dI;dR];
end

function dy=func_SIR1(t,y)
global N gamma beta1
%beta=77/(198*(3^2-1));%设接触感染者发病概率（1.19,67+10=77）
S = y(1);
I = y(2);
R = y(3);
%x=beta0*(cos(2*pi*tt)*beta1+1);
dS = -beta1*I*S/N;
dI =  beta1*I*S/N-gamma*I;
dR = gamma*I;
dy=[dS;dI;dR];
end

function dy=func_SIR2(t,y)
global N gamma beta2
%beta=77/(198*(3^2-1));%设接触感染者发病概率（1.19,67+10=77）
S = y(1);
I = y(2);
R = y(3);
%x=beta0*(cos(2*pi*tt)*beta1+1);
dS = -beta2*I*S/N;
dI =  beta2*I*S/N-gamma*I;
dR = gamma*I;
dy=[dS;dI;dR];
end

function dy=func_SIR3(t,y)
global N gamma beta3
%beta=77/(198*(3^2-1));%设接触感染者发病概率（1.19,67+10=77）
S = y(1);
I = y(2);
R = y(3);
%x=beta0*(cos(2*pi*tt)*beta1+1);
dS = -beta3*I*S/N;
dI =  beta3*I*S/N-gamma*I;
dR = gamma*I;
dy=[dS;dI;dR];
end