SIR模型简介
在典型的传染病模型中,种群(Population)内的N个个体的状态可分为如下几类:
-
易染状态S(Susceptible):即健康状态,可被感染的个体。
-
感染状态I(Infected):处于感染状态的个体还能够感染将康状态的个体。
-
移除状态R(Removed,Refractory or Recovered):也称为免疫状态和恢复状态。一个个体经历过一个完整的感染周期后,该个体就不再被感染,因此就可以不再考虑该个体。
-
传染率 β \beta β:两人接触被传染的概率,无论他是否为易感者。
-
n n n :针对于病人而言,表示了一个病人接触多少个人,可接触的人包括除自己以外的种群中的所有人
-
治愈率 γ \gamma γ:一位患者被治愈的概率
Tips:
- 初始时刻,只有少数个体处于感染状态,其他都是易染状态。
- 假设病毒的时间尺度远小于个体生命周期,从而不考虑个体的出生和自然死亡。
- 一个基本假设是完全混合(Fully mixed),也就是说一个个体与其他个体接触的机会均等。
公式
- 不考虑人口流动假设总人群数目不变: S ( t ) + I ( t ) + R ( t ) = N S(t)+I(t)+R(t)=N S(t)+I(t)+R(t)=N
- 小写字母表示某一人群数目占总人群数目的比例,大写字母表示人群数目: s ( t ) = S ( t ) N ; i ( t ) = I ( t ) N ; r ( t ) = R ( t ) N s(t)=\frac{S(t)}{N} ; i(t)=\frac{I(t)}{N};r(t)=\frac{R(t)}{N} s(t)=NS(t);i(t)=NI(t);r(t)=NR(t)
微分方程
d
S
d
t
=
−
n
β
I
(
t
)
S
(
t
)
N
\frac{dS}{dt}=-n\beta I(t)\frac{S(t)}{N}
dtdS=−nβI(t)NS(t)
d
I
d
t
=
n
β
I
(
t
)
S
(
t
)
N
−
γ
I
(
t
)
\frac{dI}{dt}=n\beta I(t)\frac{S(t)}{N}-\gamma I(t)
dtdI=nβI(t)NS(t)−γI(t)
d
R
d
t
=
γ
I
(
t
)
\frac{dR}{dt}=\gamma I(t)
dtdR=γI(t)
Tips:
- n ∗ I ( t ) n*I(t) n∗I(t)代表病人一天能接触的人总数,乘以 S ( t ) N \frac{S(t)}{N} NS(t),即易感者占人群的比例,得到患者一天接触的易感者人数,乘以感染率 β \beta β,即新增感染人数。
- 最后作图都是用的 i ( t ) , s ( t ) , r ( t ) i(t),s(t),r(t) i(t),s(t),r(t) 随 t t t 的变化的图像,所以纵坐标是一个百分比。
代码
- 感染率估计:根据解常微分方程【Matlab ODE45函数介绍】可知感染人数与时间是 y = a e b x y=ae^{bx} y=aebx关系,但是一直想不到要怎么求得感染率,最后采用武汉的确诊数据做拟合,采用matlab的拟合工具箱,最后生成拟合代码
- 移出率:找不到链接了但大致是【恢复期的倒数】这么一个说法
[num,txt,raw]=xlsread('武汉数据.xlsx');
x=num(1:22,4)';
t=1:22;
a = 50.63;
b = 0.2315;
gamma=1/14;
beta=gamma+b;
I=a*exp((beta-gamma)*t);
deta=I-x;
hold on;
plot(t,x,'r');
plot(t,I,'g');
plot(t,deta,'-o')
legend('实际确诊人数','预测确诊人数','差值')
% 自动生成 拟合代码
% %% Fit: 'untitled fit 1'.
% [xData, yData] = prepareCurveData( t, x );
%
% % Set up fittype and options.
% ft = fittype( 'exp1' );
% opts = fitoptions( 'Method', 'NonlinearLeastSquares' );
% opts.Display = 'Off';
% opts.StartPoint = [32.0107662517405 0.255589096134631];
%
% % Fit model to data.
% [fitresult, gof] = fit( xData, yData, ft, opts );
%
% % Plot fit with data.
% figure( 'Name', 'untitled fit 1' );
% h = plot( fitresult, xData, yData );
% legend( h, 'x vs. t', 'untitled fit 1', 'Location', 'NorthEast' );
% % Label axes
% xlabel t
% ylabel x
% grid on
sir.m
function y=sir(t,x)
% 感染率
a=0.3029;
% 恢复/死亡率
b=1/14;
y=zeros(3,1);
y(1)=a*x(2)*x(1)-b*x(1); %感染人数
y(2)=-a*x(2)*x(1); % 易感人数
y(3)=b*x(1); % 恢复人数
end
tspan=[0:0.5:150];% 求解区间
x01=41/10000000;
x03=3/10000000;
x02=1-x01-x03;
x0=[x01 x02 x03]; % 初值 三类人群占中人群的比例
[t,x]=ode45(@sir,tspan,x0);
hold on;
plot(t,x(:,1),'r');
plot(t,x(:,2),'g');
plot(t,x(:,3),'b');
legend('感染人数','易感人数','移出人数');
title('SIR模型图');
xlabel('天数(days)');
ylabel('总人数');