SARS传播模型.doc

学习资料收集于网络，仅供参考SARS传播模型指导老师：杜鸿飞何洪杰 刘亦珩 明 梅学习资料SARS传播模型摘要本文讨论了SARS疫情的传播规律和对经济方面的影响。首先本文对题中的早期模型进行了评价，认为其最大的优点是：能较好的描述疫情早期的发展情况，并在理论上可大致预测出疫情的爆发点以便卫生部门及早控制；最大的不足是：原模型在求解过程中参数经过多次手工调整，而且L取为一个定值，此做法主观性太强，缺乏普适性。针对上述模型的不足，本文在原模型基础上进行了改进，在非典传播的全过程中将K表示成一个函数（用Logistic函数表示），根据北京4月20日以后25个以上的数据对K进行拟合（用30个数据拟合效果较好），确定K的函数关系式，从而得到对整个过程累计病例数和日增病例变化的拟合曲线，发现它与实际情况符合得较好，而且可以再现非典传播的全过程。同时，还对K的取值进行了分析。经计算知，在相同的控制力度下，卫生部门如果提前5天采取措施，累计病例将控制在2000人以下；如果再延后5天，累计病历将至少达到3000多人，甚至可能超过4000人。最后，分析了建立一个真正能预测以及能为预防和控制提供可靠、足够信息的模型的最大困难是：因为缺乏前期数据，而不能在比较早的时候得到预测结果。本文还通过对北京市19972002年各月接待的海外旅游人数的分析并建立了时间序列模型，“预测”出2003年疫情期间本应接待的人数，对比实际接待人数，计算出在非典期间少接待的游客人数约为115万人，经济损失约1.2亿美元，约占正常情况下全年收入的33%.最后我们给当地报刊写了一篇短文，说明了建立传染病数学模型的重要性。SARS传播模型一. 问题重述SARS是21世纪第一个在世界范围内传播的传染病。SARS的爆发和蔓延给我国的经济发展和人民生活带来了很大影响，我们从中得到了许多重要的经验和教训，认识到定量的研究传染病的传播规律、为预测和控制传染病蔓延创造条件的重要性。所以本文首先评价了一个已有的早期模型的合理性和实用性，然后在此基础上建立了一个更优的模型并给出了分析，最后建立了SARS对经济某个方面影响的模型。二.问题分析（一）通过对早期模型和实际情况的分析，我们认为影响SARS传播因素众多，大致可分为时域因素和地域因素。列举如下：(1)时域因素a媒体宣传：初期疫情较轻，媒体宣传强度很弱，导致民众的自我保护意识不足，容易感染；后期疫情较重，媒体宣传强度很大，民众的自我保护意识大大加强。b政府干预：初期疫情较轻，政府介入不足，后期疫情较重，政府加强干预（如：强行隔离，公共场所消毒等行为）。c认识程度：当一种新的传染病出现时，初期由于人们的认识程度不足，无法采取有效的预防和治疗措施，但随着研究的深入，认识程度会越来越高。(2)地域因素a经济水平和医学水平：经济水平和医学水平高的地区的疫情控制情况会明显比水平低的地方好。b人口密度和人口流动：人口密度和人口流动大的城市若爆发传染病，疫情程度会比人口密度和人口流动小的城市大。c气候：SARS适合在春秋两季传播，且各城市的气候会疫情程度。（二）我们认为在SARS疫情期间考察的人群大致可分为三类：健康人群，感染人群，治愈人群。而感染人群又可分为非传染源和传染源两类。（三）一个较好的传染病传播模型因该具有如下功能： a能较好的描述疫情的大致走势。 b能较精确的给出关键时间（初期爆发时刻；中期稳定时刻；高峰期；0病例增长的时刻），以便政府和卫生部门针对不同作出及时而正确的措施。 c能给出描述疫情的指标，以便政府和卫生部门决定其各项工作的力度。三基本假设1题中所给的数据真实可信。2假定疫情爆发后政府一定会采取措施。3假定北京市医院及医务人员足够多。四变量说明N(t) 累计病例数K 平均每病人每天可传染人数L 平均每个病人可以直接感染他人的天数五早期模型评价（一）早期模型重述假定初始时刻的病例数为N0，平均每病人每天可传染K个人，平均每个病人可以直接感染他人的时间为L天。则在L天之内，累计病例数N(t)随时间t（单位天）的关系是： （1-1）如果不考虑传染期的限制，病例数将按指数规律增长，考虑传染期的限制后，则采用半模拟循环计算的办法，把达到L天的病例从可以引发直接传染的基数中去掉。然后假定从开始至高峰期间均采用同样的K值（从拟合这一阶段的数据得出），到达高峰期后，在10天的范围内逐步调整K值到比较小，然后保持不变，拟合其后在控制阶段的全部数据（认为社会在短期剧烈调整后，进入对疫情控制较好的常态）。（二）早期模型的合理性和实用性的评价A早期模型的优点1模型（1-1）实际上是微分方程在（0L）区间内的特解1。其中N(t)表示t时刻的累计病例数，则(N(t)-N(t-L)表示传染源数量，为病例总数减去失去传染能力的病例数。2参数K和L是描述SARS传播的两个重要参数，并且具有实际的意义：L可理解为平均每个病人在被发现前后可以造成造成直接传染的期限，在此期限后失去传染作用，可能原因是被隔离、病愈或死去等等。K表示某种社会条件下平均每病人每天传播的人数（但并非文中所述的一个病人的感染他人的平均概率）。3通过对模型的分析，可得到一预测疫情发展的参数RKL，R表示平均每个病人在其传染期内感染的总人数。若R1，说明社会上现存传染源人数在上升，疫情将因失控而爆发；若R1，说明社会上现存传染源人数在下降，疫情将得到控制；若R1，说明社会上现存传染源人数不变，疫情将持续下去2。4从拟合的图形来看，此模型对各城市早期的SARS疫情描述的较好，具有一定的通用性。其实际意义就是此模型可大致预测出疫情的爆发点和发展趋势。因为通过对早期数据的拟合确定参数，得出形如（1-1）的一个指数形式的模型，而指数函数的曲线初期增长较慢，后期增长急速，必可大致找到一个“转折点”，而“转折点”所对应的时间便是预测的疫情爆发点。这一数据对于卫生部门十分重要，因为控制疫情的最好时间是在疫情爆发之前。B早期模型的不足之处1首先模型并未给出拟合程度的参数，而当我们试图通过计算得到该模型的拟合程度参数时发现无法进行。原因是原模型求解过程的中间阶段参数K多次手工调整，而且模型中并未给出调整的标准和相关理论，所以我们很难重复该求解过程。由此我们得出结论：该模型的参数取值主观性太强，此作法给阅读者运用并改进模型带来了极大的困难，所以此模型的普适性较差。2在数据不足的情况下因无法进行手工调整，所以该模型用香港后期拟合的K值去预测北京后期疫情的发展趋势。但如问题分析所述，地域因素会造成不同地区的K值不同（如人口密度和人口流动大的城市若爆发传染病，初期的K值会比人口密度和人口流动小的城市大，等等），而很难找到地域因素几乎相同的两城市。所以此作法可能导致预测结果相差较大。图1为用此方法预测的北京后期疫情情况与实际情况的对比图 图1从图中可以看出，预测值与真实值偏差越来越大。对该模型的评价：该模型具有较好的实际意义，能比较合理的反映非典的传播情况和发展趋势，模型中的参数设置也较科学，能较好地反映非典传播的影响因素；但模型的求解过程主观性太强，很难重现，而且参数的取值也很主观，没有取值的标准和理论。综上所述，我们得到结论：该模型较合理，但不实用。六SARS传播模型的建立1模型的建立与参数讨论：假设N(t)为随时间变化的累积SARS病人总数；K(t)为某一天平均每病人传染他人的人数（/天），是时间的函数；L为平均每个病人可传染他人的传染期限（天）；则为单位时间（天）内增加的发病的人数，表示第t天时具有传染能力的人数，K(t)(N(t)-N(t-L)就表示t天时增加的被感染的人数。由上分析得到N(t)的微分方程如下： （1）同时，随着时间的变化，由于外界因素的改变，人均日传染数K(t)也会变化：在疫情初发期，由于人们对SARS并没有什么认识，更不知道其严重性，所以即使有患病者，该患病者的活动范围也比较大，可能传染的人数也比较多，故人均日传染数K(t)就比较大并且在一定时间内保持基本稳定；当患者越来越多，疫情越来越严重时，SARS就会受到社会的普遍关注，政府部门也会立即采取强制控制手段来限制疫情的发展，民众自我防范意识在媒体宣传等作用下加强，患病和疑似患病者活动范围受到严格控制，从而人均日传染数K(t)开始快速下降；当人均日传染数K(t)下降到一个较小值之后，由于传染的可能性仍存在，政府部门的控制能力也有一定限度，K(t)的下降速度明显变缓；最后随着累积病例数的稳定，K(t)缓慢降至0。通过以上分析，可以得到如下的类Logstic微分方程：， （2）其中M为K的一个上界，r为衰减系数，它表示K的变化率与K(M-K)成反比，一方面，K的变化率与K本身的大小有关：当K值很大时，疫情较严重，无论政府干预行为还是民众自我保护行为都很强，所以K下降很快；当K很小时，情况相反。另一方面，K的变化率与的(M-K)大小有关：当(M-K)很小时，K的下降“空间”很大，对政府干预和民众保护行为很“敏感”，下降的很快；当(M-K)很大时，K的下降“空间”很小，这时情况相反。图2就反映了我们上面分析的K的变化趋势：图2 人均日传染数K随时间（3月1号起）的变化图参数L可理解为平均每个病人可以造成直接传染的期限，在此期限后他失去传染作用，可能的原因是被严格隔离、病愈不再传染或死去等等。与题目附件一中的讨论一样，我们认为L与医疗机构有关，取L=20这个具有一定统计意义的值。综上，我们得到了SARS传染的微分方程模型如下：，且，L=20（其中“t=1”指出现第一例病人的时间。）2模型的求解：A模型推导求解方程（2），得到人均日传染数K(t)的函数关系式如下： ， 其中且其中M为K的一个上界，衰减系数r0，并且由前面的讨论，知当时应有，故应有参数cdd mindd=dd; %记录最小的D和相应的M，c，r markM=M;markc=c;markr=r;endendmarkM,markc,markr %找到M，c，r%程序2 代入M计算拟合与预测的N(t)值M=markM;c=markc;r=markr; n0=9;for t=1:160 k(t)=M/(1-c*exp(M*r*t); %由表达式求出K的所有值endfigure(1);plot(k)%求拟合与预测的累计病例数N(t)之值n(1)=1;for t=2:21 n(t)=n(t-1)*(1+k(t);endfor t=22:115 n(t)=n(t-1)+(n(t-1)-n(t-20-1)*k(t);end%由上面的N(t)计算日增量dN(t)for t=1:114 dn(t)=n(t+1)-n(t);endfigure(2); %画图plot(n,ro)hold onra(51:115)=a(1:65);plot(51:115,ra(51:115),b*)hold offfigure(3);plot(dn,)附件2：问题3的程序%程序1 用傅立叶级数进行拟合%问题3的数据y0=9.4 11.3 16.8 19.8 20.3 18.8 20.9 24.9 24.7 24.3 19.4 18.6 9.6 11.7 15.8 19.9 19.5 17.8 17.8 23.3 21.4 24.5 20.1 15.9 10.1 12.9 17.7 21.0 21.0 20.4 21.9 25.8 29.3 29.8 23.6 16.5 11.4 26.0 19.6 25.9 27.6 24.3 23.0 27.8 27.3 28.5 32.8 18.5 11.5 26.4 20.4 26.1 28.9 28.0 25.2 30.8 28.7 28.1 22.2 20.7 13.7 29.7 23.1 28.9 29.0 27.4 26.0 32.2 31.4 32.6 29.2 22.9;%进行傅立叶级数拟合w=2*pi/12;for t=1:72 ac(t)=cos(2*w*t);ac1(t)=cos(w*t); bc(t)=sin(2*w*t);bc1(t)=sin(w*t);endx=ac;ac1;bc;bc1;1:72;ones(1,72);b=xy0;%拟合得到的表达式的系数%b=-2.2503 -2.7344 -2.6120 -3.5255 0.1545 16.7552%计算出拟合出的傅立叶级数的值t=1:72;y=b(1)*cos(2*w*t)+b(2)*cos(w*t)+b(3)*sin(2*w*t)+b(4)*sin(w*t)+b(5)*t+b(6)*ones(1,72);%计算出拟合出的傅立叶级数的值并计算出未来8个月的值tt=1:80;yy1=b(1)*cos(2*w*tt)+b(2)*cos(w*tt)+b(3)*sin(2*w*tt)+b(4)*sin(w*tt)+b(5)*tt+b(6)*ones(1,80);yy=yy1(73:80);%程序2 用时间序列拟合%计算真实值与傅立叶级数计算出的值的差z=y0-y;%计算模型的参数的值r0=sum(z.2)/72;t1=0;for i=1:71 t1=t1+z(i)*z(i+1);endr1=t1/72;t2=0;for i=1:70 t2=t2+z(i)*z(i+2);endr2=t2/72;p1=r1/r0;p2=r2/r0;phi1=p1*(1-p2)/(1-p12);phi2=p2*(1-p1)/(1-p12);%计算方差，但此处数据较少，方差值可能不准确segma=sqrt(r0*(1-phi1*p1-phi2*p2);%用题目中的数据检验预测的