传染病SIR模型

SIR模型简介

在典型的传染病模型中，种群（Population）内的N个个体的状态可分为如下几类：

• 易染状态S(Susceptible)：即健康状态，可被感染的个体。

• 感染状态I(Infected)：处于感染状态的个体还能够感染将康状态的个体。

• 移除状态R(Removed,Refractory or Recovered)：也称为免疫状态和恢复状态。一个个体经历过一个完整的感染周期后，该个体就不再被感染，因此就可以不再考虑该个体。

• 传染率$\beta$：两人接触被传染的概率，无论他是否为易感者。

• $n$ ：针对于病人而言，表示了一个病人接触多少个人，可接触的人包括除自己以外的种群中的所有人

• 治愈率$\gamma$：一位患者被治愈的概率

Tips:

1. 初始时刻，只有少数个体处于感染状态，其他都是易染状态。
2. 假设病毒的时间尺度远小于个体生命周期，从而不考虑个体的出生和自然死亡。
3. 一个基本假设是完全混合（Fully mixed）,也就是说一个个体与其他个体接触的机会均等。

公式

• 不考虑人口流动假设总人群数目不变：$$S(t)+I(t)+R(t)=N$$
• 小写字母表示某一人群数目占总人群数目的比例，大写字母表示人群数目：$$s(t)=\frac{S(t)}{N};i(t)=\frac{I(t)}{N}；r(t)=\frac{R(t)}{N}$$

微分方程

$$\frac{dS}{dt}=-n\beta I(t)\frac{S(t)}{N}$$
$$\frac{dI}{dt}=n\beta I(t)\frac{S(t)}{N}-\gamma I(t)$$
$$\frac{dR}{dt}=\gamma I(t)$$

Tips：

• $n\ast I(t)$代表病人一天能接触的人总数，乘以$\frac{S(t)}{N}$，即易感者占人群的比例，得到患者一天接触的易感者人数，乘以感染率$\beta$，即新增感染人数。
• 最后作图都是用的$i(t)，s(t),r(t)$ 随 $t$ 的变化的图像，所以纵坐标是一个百分比。

代码

• 感染率估计：根据解常微分方程【Matlab ODE45函数介绍】可知感染人数与时间是$y=a{e}^{bx}$关系，但是一直想不到要怎么求得感染率，最后采用武汉的确诊数据做拟合，采用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('总人数');