%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