传染病SIR模型及蒙特卡洛方法 |
您所在的位置:网站首页 › 传染病公式 › 传染病SIR模型及蒙特卡洛方法 |
传染病模型的Matlab实现
一、符号说明二、基本公式三、基本数据处理四、rβ参数估计及模型简化五、SIR模型MATLAB实现及分析六、蒙特卡洛方法的实现及分析七、模型改进八、文章说明
一、符号说明
以SIR模型为例,这里对SIR重新定义: S:Susceptible 易感人群(数量),与感染者密切接触易受到其感染。I:Infective 感染人群(数量),已患上传染病并有可能传染给易感人群。由于在我国内,感染者一旦出现症状便被隔离,几乎丧失传染的可能性,故被隔离患者归入移除人员,即不再参与传染过程的人员。R:Removed 移除人群(数量),患病后被治愈或死亡人群数量,这部分人不再参与感染和被感染过程。由于国家的隔离政策,被隔离人群也被计算在失去传染能力的移除人群内。N: 人口总数,包括易感人群,感染人群及移除人群。r: 每个患者每天内接触在易感者人数。β: 易感者与感染者接触时的传染率。γ: 感染人群失去传播能力的速率,其中因隔离而丧失传播能力的部分占比应远远大于自行恢复及死亡的部分。注:(rβ)可以反映出每一个患者的传染能力。反映了疾病的传染强度,其值越大说明感染人群与易感染人群接触时传染的可能性越大。r可以通过统计方法进行计算,而为了追求实际意义,我们后文主要将(rβ)看作一个整体。 . . 二、基本公式转换关系的简单示意图如下: 我们为了简化方程,假设人口基数足够大,即在某一较短时间S/N的改变几乎可以忽略不计,举疫情传播早期为例,在疫情刚开始I/N≈0,在短期内S/N≈(N-I)/N≈1,则式(2)可简化为: 数据来源:839Studio. Novel-Coronavirus-Updates [DB/OL].GitHub,10:30 May. 8 GMT+8 我们通过将数据处理得到武汉市1月18日至2月1日各类数据,选择这期间的数据的原因是,1月18日接近疫情初期,而2月1日为武汉封城日期。 . 四、rβ参数估计及模型简化利用MATLAB非线性拟合模块确定参数rβ粗略值,编写MATLAB程序如下: load WuhanData.mat %导入1月18日至2月1日武汉各类人群数据 %数据格式: %====================================== %月份 日期 累计确诊 累计出院 累计死亡 % 1 18 45 15 2 % 1 19 62 19 2 [m,n]=size(WuhanData); WuhanData(:,3)=WuhanData(:,3)-WuhanData(:,4)-WuhanData(:,5); %I=e^((rβ-1/14)t) ft=fittype('n*exp((lmd-1/14).*t)','coefficients',{'lmd','n'},... 'independent','t','dependent','I'); cf=fit((1:m)',WuhanData(:,3),ft,'StartPoint',[0 41]); plot(cf,(1:m)',WuhanData(:,3)) disp(cf)解得: General model: cf(t) = n*exp((lmd-1/14).*t) Coefficients (with 95% confidence bounds): lmd = 0.3123 (0.273, 0.3516) n = 104.5 (52.87, 156.1) 我们将使用MATLAB进行数值模拟,来向大家展示各个参数对病毒传播的影响。 为了简化模型,忽略地区人数对实验结果造成的影响,仅对模型的本质进行模拟,并将各类人群以比例的方式表示,即将原本以个人为单位的离散化模型变得连续化,我们可假设我们所探究区域的人口总数在短期内基本保持不变,即:
编写SIR模型代码如下: [t,sir]=ode15s('SIR',[0 250],[1-1e-6;1e-6;0]); hold on plot(t,sir(:,1),'linewidth',1.5) plot(t,sir(:,2),'linewidth',1.5) plot(t,sir(:,3),'linewidth',1.5) [~,peakTpos]=max(sir(:,2)); scatter(t(peakTpos),sir(peakTpos,2),20,'k','filled') legend('S:易感者比例','I :感染者比例','R:移除者比例','感染者峰值坐标') title('SIR模型:S(0)=1-10^{-6} I(0)=10^{-6} r\beta=0.3123 \gamma=1/14') xlabel('t:时间(天)') ylabel('各类人群所占比例') function newrate=SIR(t,rate) %rate: %================== %[s(t) i(t) r(t)] rb=0.3123; g=1/14; newrate=[-rb*rate(1)*rate(2);... rb*rate(1)*rate(2)-g*rate(2);... g*rate(2)]; end若取初始感染者占比为
1
0
−
6
10^{-6}
10−6,即假设该市有1000万人口,在疫情初始阶段约有10人患病,求解常微分方程数值解结果如下: 我们假设(rβ)中r=5,即每人每天平均接触五个人,我们如果进行自我隔离,将r分别降至2或1,得到结果如下: 假设我们加强了检测力度、拥有了更加先进的对该病毒的检测方法、或对该病毒的医疗有了突破性进展,我们将γ提高至1/7,则的得到如下结果: 通过各种手段降低β的情形与降低r类似,由于该连续模型下,(rβ)被当作一个整体看待,很难将其分别进行分析,故我们在此引入蒙特卡洛方法: 在一定区域内,将每个人看作独立的个体,用圆点来表示,每隔单位时间每个个体都会向一定方向位移一段距离并且当易感者和感染而未被隔离者相距较近时,易感者可能会向感染而未被隔离者转换,我们使用灰色圆点表示易感者,粉色圆点表示感染而未被隔离者,红色圆点表示移除或被隔离者,红色圆点已不具备传染能力。编写蒙特卡洛方法代码如下: function MonteCarlo figure('units','pixels',... 'position',[350 100 800 500],... 'Numbertitle','off',... 'MenuBar','none',... 'resize','off',... 'name','MonteCarlo',... 'color',[0.95 0.95 0.95]); axes('position',[0 0 1 1]); axis([0 800 0 500]); set(gca,'color',[0 0 0.0078]); set(gca,'xtick',[],'ytick',[],'xcolor','w','ycolor','w') %========================================================================== Limit=[0 800 0 500]; beta=0.3; gama=0.01; InfectDis=1; Dis=2; Num.N=500; Num.S=Num.N-5; Num.I=5; Num.R=0; Pos.N=rand(Num.N,2).*[Limit(2),Limit(4)]; Pos.S=Pos.N(1:Num.S,:); Pos.I=Pos.N(Num.S+1:Num.S+Num.R,:); Pos.R=Pos.N(Num.S+Num.R+1:end,:); Ncollection=Num.N; Scollection=Num.S; Icollection=Num.I; Rcollection=Num.R; %========================================================================== hold on; DrawSHdl=scatter(Pos.S(:,1),Pos.S(:,2),50,'filled','CData',[0.3,0.3,0.3]); DrawIHdl=scatter(Pos.I(:,1),Pos.I(:,2),50,'filled','CData',[1,0.5,0.5]); DrawRHdl=scatter(Pos.R(:,1),Pos.R(:,2),50,'filled','CData',[1,0.1,0.1]); TxtNHdl=text(10,480,['人口总数:',num2str(Num.N)],'color',[1 1 1],'fontsize',12,'FontWeight','bold'); TxtIHdl=text(10,480-25,['携带者:',num2str(Num.I)],'color',[1 1 1],'fontsize',12,'FontWeight','bold'); TxtRHdl=text(10,480-50,['隔离或具抗性者:',num2str(Num.R)],'color',[1 1 1],'fontsize',12,'FontWeight','bold'); fps = 200; game = timer('ExecutionMode', 'FixedRate', 'Period',1/fps, 'TimerFcn', @Spread); set(gcf,'tag','co','CloseRequestFcn',@clo); start(game) function clo(~,~) stop(game) save SIR.mat Ncollection Scollection Icollection Rcollection delete(findobj('tag','co')); clf close end function Spread(~,~) Theta=rand(Num.N,1).*2.*pi.*Dis; MovDir=[cos(Theta),sin(Theta)]; Pos.N=Pos.N+MovDir; Pos.S=Pos.N(1:Num.S,:); Pos.I=Pos.N(Num.S+1:Num.S+Num.I,:); Pos.R=Pos.N(Num.S+Num.I+1:end,:); DisMat=GetDis(Pos.S,Pos.I); S2I=any(DisMat |
今日新闻 |
推荐新闻 |
CopyRight 2018-2019 办公设备维修网 版权所有 豫ICP备15022753号-3 |