%f(1)的起点:
f(1:7)=[1,2,3,4,5,6,7]; %初始化自由带菌者数。 bb=0.9;k=3; for i=1:N
a(i+6)=f(i)+a(i+6); func1(k*f(i),i+6); func2(k*f(i),i+7); end
bb=.4;k=1; for i=N:180
a(i+6)=f(i)+a(i+6); func1(k*f(i),i+6); func2(k*f(i),i+7); end
P=P+a; end P=P/M;
nn=length(find(P>0.1)); figure;
AA=[60,70,80,100,143,106,105,81,103,111,126,85,148,93,113,83,105,62,94,63,89,87,41,50,38,...
39,43,23,18,17,15,14,3,7,0,12,9,25,9,5,8,2,3,3,1,1,0,0,0,0,0,1,-1,0,0,1,0,-1,0,0,-1,0,0,... 0,0,0,0,0];
plot(37:36+length(AA),AA,'r:');hold; plot(1:nn,P(1:nn),'b--');hold legend('实际','拟和') title('改进后的模拟')
function func2(q,p) global a f bb
%func1---指q人群被感染后当天(感染源就在当天发现)就开始被寻找。(成为隔离的人群)
%b是遗漏率。0
%1-b(1),b(1)*(1-b(2)),b(1)*b(2)*(1-b(3))
bbb=bb^(1/3); ?b为每天的遗漏率,bb为总遗漏率。b是模拟值。 if q<.5 else
n=floor(.5+normrnd(5,2)); if n<1;n=1;end for i=1:3
19
b(i)=normrnd(bbb,.05); if b(i)<0 b(i)=0; elseif b(i)>1 b(i)=1; end end
if n>2
a(p+n)=(1-b(1)*b(2)*b(3))*q; %一部分染病者在被隔离时发病
f(p+n)=b(1)*b(2)*b(3)*q; %另一部分未被隔离,成为自由带菌者
elseif n==2
a(p+n)=(1-b(1)*b(2)*b(3))*q; f(p+n)=b(1)*b(2)*b(3)*q;
func2(b(1)*b(2)*(1-b(3))*q,p+2); %第二天虽被找到,但已造成当天的影响。
elseif n==1
a(p+n)=(1-b(1)*b(2))*q; %已在隔离区中(第零,一天找到的人数和)
a(p+n+1)=b(1)*b(1)*q; %迟一天才找到的(第一天之前(发病之前)未被隔离)。
func1(b(1)*b(2)*q,p+2);%他从发病到入院,有两天的自由感染期。(第二天才被找到)
func2(b(1)*b(2)*q,p+3);
func2(b(1)*(1-b(2))*q,p+2); %第一天虽被找到,但已造成当天的影响。 end end
function func1(q,p) global a f bb
%q是接触患者后受到感染的人数,p是他们被感染的日期。
%func1---指q人群被感染后一天(感染源才被发现)才开始被寻找,(成为隔离的人群).
%b是遗漏率。0
%1-b(1),b(1)*(1-b(2)),b(1)*b(2)*(1-b(3))
bbb=bb^(1/3); ?b为每天的遗漏率,bb为总遗漏率。b是模拟值。 if q<.5 else
n=floor(.5+normrnd(5,2)); if n<1;n=1;end
20
for i=1:3
b(i)=normrnd(bbb,.05); if b(i)<0 b(i)=0; elseif b(i)>1 b(i)=1; end end if n>3
a(p+n)=(1-b(1)*b(2)*b(3))*q; f(p+n)=b(1)*b(2)*b(3)*q; elseif n==3
a(p+n)=(1-b(1)*b(2)*b(3))*q; f(p+n)=b(1)*b(2)*b(3)*q;
func2(b(1)*b(2)*(1-b(3))*q,p+3); %第三天虽被找到,但已造成一天的影响。
elseif n==2
a(p+n)=(1-b(1)*b(2))*q; %已在隔离区中(第一,二天找到的人数和)
a(p+n+1)=b(1)*b(1)*q; %迟一天才找到的(第二天之前(发病之前)未被隔离)。
func1(b(1)*b(2)*q,p+2); %他从发病到入院,有两天的自由感染期。
func2(b(1)*b(2)*q,p+3);
func2(b(1)*(1-b(2))*q,p+2); %第二天虽被找到,但已造成一天的影响。
elseif n==1
a(p+n)=(1-b(1))*q; %分析同上。 a(p+n+1)=b(1)*q; func1(b(1)*q,p+1); func2(b(1)*q,p+2);
func2((1-b(1))*q,p+1); end end
附件3 SARS减少的每日旅游人数 Clear;
Y=zeros(1,160);
for i=1:160;%以-0.00976159133019为斜率 0.00684210526316为切距 Y(i)=-0.00976159133019*i+0.684210526316; end
plot(1:160,Y);%8月15日后SARS对每日旅游人数的减少量,10月25后影响消
21
失。
a=sum(Y(1,1:15))%8月15到8月底 b=sum(Y(1,16:45))%9月受减小量 c=sum(Y(1,46:70))月受减小量。
22