传染病模型：
常见的传染病模型按照具体的传染病的特点可分为 SI、SIS、SIR、SIRS、SEIR 模型。其中“S”“E”“I”“R”的现实含义如下：
S (Susceptible)，易感者，指缺乏免疫能力健康人，与感染者接触后容易受到感染；
E (Exposed)，暴露者 ，指接触过感染者但暂无传染性的人，可用于存在潜伏期的传染病；
I (Infectious)，患病者，指有传染性的病人，可以传播给 S，将其变为 E 或 I ；
R (Recovered)，康复者，指病愈后具有免疫力的人，如是终身免疫性传染病，则不可被重新变为 S 、E 或 I ，如果免疫期有限，就可以重新变为 S 类，进而被感染。
SIS传染病模型假设：
1.在疾病传播期内，研究区域内总人数N不变,既不考虑生死，也不考虑迁移。将该区域内的人群分为：易感染者、感病者和病愈免疫的移出者三类，以下简称易感染者，病人和移出者。在t时刻将这三类人群的数量分别记为S(t)、I(t)和R(t)，t时刻这三类人在总人数中所占的比例分别记为s(t)、i(t)和r(t)。
2.每个病人每天有效接触的平均人数是常数λ，称为日接触率。当病人与健康者有效接触时,会使健康者受感染变为病人。
3.每天被治愈的病人占病人总数的比例为常数μ，称为日治愈率。病人被治愈后仍有可能被感染为病人，故假设1/μ是该疾病的平均传染期。
4.在初始时刻，只有少数个体处于感染状态，其他都是易感染状态。
5.假设疾病的时间尺度远小于个体生命周期，从而不考虑个体的出生和自然死亡对数据统计造成的影响。
6.完全混合：每一个个体与其他个体接触的机会均等。
根据以上假设可得：每个病人每天可使λ·s(t)个健康者变为病人，每天有μ·Ni个病人被治愈。
σ=λ/μ
由该式中，接触数σ是病人在平均传染期内有效接触的人数。
SIS传染病模型构建：
由假设1可得，整个群体由易感染者，感染者，移出者构成，三者之间的关系如下：
根据假设1可得：
由上式可得，关于S(t),I(t),R(t)的微分方程：
进而得到关于s(t),i(t),r(t)的微分方程：
将数据带入上述公式，求解微分方程即可得到各个变量随时间变化的情况。
matlab代码：
function [t,S,I] = Program_7_7(N,n,tau,gamma,MaxTime,Type)
% It is an SIS disease spread through a network. Allowed
% network types are 'Random','Spatial','Lattice' and 'SmallWorld'
%
% We assume N individuals, each with an averge of n contacts.
% In this model we define an individual by their status flag:
% Status=1 => Susceptible
% Status=2 => Infectious
% Status=0 => Recovered
% Sets up default parameters if necessary.
if nargin == 0
N=100;
n=4;
tau=1;
gamma=0.1;
MaxTime=1000;
Type='Random';
end
% Checks all the parameters are valid
CheckGreater(N,0,'Number of individuals N');
CheckGreater(n,0,'Number of neighbours n');
CheckGreater(tau,0,'tau');
CheckGreater(gamma,0,'gamma');
CheckGreater(MaxTime,0,'MaxTime');
%Initialise the Network
% (X,Y) is location, G is the network graph matrix
% this means we use S and I for the number of susceptibles and infecteds
[X,Y,G,N]=Create_Network(N,n,Type);
Status=1+0*X; Status(1)=2;
Rate=0*X; Rate(1)=gamma; Rate(find(G(:,1)))=tau;
t=0; i=1; S=N-1; I=1;
% The main iteration
subplot(2,1,1);
[j,k,s]=find(G);
plot([X(j) X(k)]',[Y(j) Y(k)]','-k');
hold on
Col=[0.7 0.7 0.7; 0 1 0; 1 0 0];
for k=1:N
H(k)=plot(X(k),Y(k),'ok','MarkerFaceColor','g');
end
set(H(1),'MarkerFaceColor','r');
hold off;
drawnow;
while (t<MaxTime & I(end)>0)
[step,Rate,Status,e]=Iterate(Rate,Status,G,tau,gamma);
i=i+1;
t(i)=t(i-1)+step;
S(i)=length(find(Status==1));
I(i)=length(find(Status==2));
set(H(e),'MarkerFaceColor',Col(Status(e)+1,:));
% subplot(2,1,1);
% [j,k,s]=find(G);
% plot([X(j) X(k)]',[Y(j) Y(k)]','-k');
% hold on
% s=find(Status==0); plot(X(s),Y(s),'ok','MarkerFaceColor',[0.7 0.7 0.7],'MarkerSize',8);
% s=find(Status==1); plot(X(s),Y(s),'ok','MarkerFaceColor','g','MarkerSize',8);
% s=find(Status==2); plot(X(s),Y(s),'ok','MarkerFaceColor','r','MarkerSize',10);
% hold off;
subplot(4,1,3);
plot(t,S,'-g');
ylabel 'Number of Susceptibles'
subplot(4,1,4);
plot(t,I,'-r');
ylabel 'Number of Infecteds'
xlabel 'Time'
drawnow;
[];
end
% Create the Network
function [X,Y,G,N]=Create_Network(N,n,Type);
if n > (N-1)
error('Impossible to have an average of %d contacts with a population size of %d',n,N);
end
G=sparse(1,1,0,N,N);
X=rand(N,1); Y=rand(N,1);
switch Type
case {'Random','random'}
contacts=0;
while(contacts<n*N)
i=ceil(rand(1,1)*N); j=ceil(rand(1,1)*N);
if i~=j & G(i,j)==0
contacts=contacts+2;
G(i,j)=1; G(j,1)=1;
end
end
case {'Spatial','spatial'}
D=(X*ones(1,N) - ones(N,1)*X').^2 + (Y*ones(1,N) - ones(N,1)*Y').^2;
Prob=exp(-D*5)-rand(N,N); Prob=triu(Prob,1)-1e100*tril(ones(N,N),0);
[y i]=sort(reshape(Prob,N*N,1));
p=i(end+[(1-n*N/2):1:0]);
i=1+mod(p-1,N); j=1+floor((p-1)/N);
G=sparse([i; j],[j; i],1,N,N);
case {'Lattice','lattice'}
if mod(sqrt(N),1)~=0
warning('N=%d is not a square number, rounding to %d',N,round(sqrt(N)).^2);
N=round(sqrt(N)).^2;
end
[X,Y]=meshgrid([0:(sqrt(N)-1)]/(sqrt(N)-1));
X=reshape(X,N,1); Y=reshape(Y,N,1);
D=(X*ones(1,N) - ones(N,1)*X').^2 + (Y*ones(1,N) - ones(N,1)*Y').^2;
Prob=triu(D,1)+1e100*tril(ones(N,N),0);
[y i]=sort(reshape(Prob,N*N,1));
p=i(1:(n*N/2 - 2*sqrt(N)));
i=1+mod(p-1,N); j=1+floor((p-1)/N);
G=sparse([i; j],[j; i],1,N,N);
case {'SmallWorld','Smallworld','smallworld'}
if mod(sqrt(N),1)~=0
warning('N=%d is not a square number, rounding to %d',N,round(sqrt(N)).^2);
N=round(sqrt(N)).^2;
end
[X,Y]=meshgrid([0:(sqrt(N)-1)]/(sqrt(N)-1));
X=reshape(X,N,1); Y=reshape(Y,N,1);
D=(X*ones(1,N) - ones(N,1)*X').^2 + (Y*ones(1,N) - ones(N,1)*Y').^2;
Prob=triu(D,1)+1e100*tril(ones(N,N),0);
[y i]=sort(reshape(Prob,N*N,1));
p=i(1:(n*N/2 - 2*sqrt(N)));
i=1+mod(p-1,N); j=1+floor((p-1)/N);
G=sparse([i; j],[j; i],1,N,N);
% Already made a lattice, now add R random connections
R=10;
for k=1:R
i=1; j=1;
while (i==j | G(i,j)==1)
i=ceil(rand(1,1)*N); j=ceil(rand(1,1)*N);
end
G(i,j)=1; G(j,1)=1;
end
otherwise
error('Do not recognise network type %s',Type);
end
%Do the Up-Dating.
function [step,Rate,Status,Event]=Iterate(Rate,Status,G,tau,gamma);
Sum=sum(Rate); Cum=cumsum(Rate);
Event=min(find(Cum>rand(1,1)*Sum));
Status(Event)=mod(Status(Event)+1,3);
contacts=find(G(:,Event) & Status==1);
switch Status(Event)
case 1
case 2
Rate(Event)=gamma;
Rate(contacts)=Rate(contacts)+tau;
G(Event,:)=0;
case 0
% For SIR type dynamics we require the following 2 lines
Rate(Event)=0;
Rate(contacts)=Rate(contacts)-tau;
% For SIS type dynamics we require the following 3 lines
%Status(Event)=1;
%Rate(Event)=tau*length(find(G(:,Event) & Status==2));
%Rate(contacts)=Rate(contacts)-tau;
end
Rate=Rate.*sign(Status);
step=-log(rand(1,1))/Sum;
% Does a simple check on the value
function []=CheckGreaterOrEqual(Parameter, Value, str)
m=find(Parameter<Value);
if length(m)>0
error('Parameter %s(%g) (=%g) is less than %g',str,m(1),Parameter(m(1)),Value);
end
function []=CheckGreater(Parameter, Value, str)
m=find(Parameter<=Value);
if length(m)>0
error('Parameter %s(%g) (=%g) is less than %g',str,m(1),Parameter(m(1)),Value);
end
结果
1022
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
