传染病模型

最新推荐文章于 2023-06-03 14:36:14 发布

无人机_Fly

最新推荐文章于 2023-06-03 14:36:14 发布

阅读量1k

点赞数 2

CC 4.0 BY-SA版权

文章标签：机器学习算法 python

本文链接：https://blog.youkuaiyun.com/qq_60337394/article/details/126131578

博客内容探讨了不同的传染病模型，包括SI、SIS、SIR和更复杂的SIRS、SEIR模型，用于模拟病毒如何在人口中传播。通过调整参数如感染率、治愈率和复感率，模拟不同情况下疾病的传播动态。此外，还提供了一个具体的MATLAB代码示例，模拟了用户网络中病毒传播的过程，考虑了用户的预免疫、病毒变异和用户对病毒的反应时间等因素，展示了不同设置下易感者、潜伏者、患病者和康复者数量的变化。

摘要生成于 C知道，由 DeepSeek-R1 满血版支持，前往体验 >

S（易感者）：缺乏免疫力的健康人，与感染者接触容易感染；
E（暴露者/潜伏期）：接触过感染者或者感染后恢复健康的人，此时并未变成患病者，但不具有传染性；
I （患病者）：需要接受治疗且具有传染性；
R（康复者）：病愈后具有免疫力的人，但不代表没有变成S的可能性。

大部分的传染病模型都是关于时间t的函数

SI模型

不会反复发作的疾病
总人数为N,N=I+S,记i和s为易感者S和患病者I占总数N得比例，S=N $⋅\cdot$ s.
$λ\lambda$ 是单个患病者（携带病毒可以传染给他人的人）每天能接触到的易感者比例。
患病者每天增加的比例：
$didt=i×λ⋅s=i×(1−i)λ\frac{di}{dt}=i\times \lambda \cdot s = i \times (1-i)\lambda$
患病者每天增加的数量：
$dIdt=I×λ⋅S/N=I×(1−I)λ/N\frac{dI}{dt}=I\times \lambda \cdot S/N= I\times (1-I)\lambda /N$
求解得原函数:
$I(t)=N1+(N/I0−1)⋅e−λtI(t)=\frac{N}{1+(N/I_0 - 1)\cdot e^{-\lambda t}}$


clc;
close all;
clear all;
I=10;
N=10000;
S=N-I;
lemda=0.1;
t=1:365;
for i=1:(size(t,2)-1)
    I(1+i)=I(i)+I(i)*(N-I(i))*lemda/N;
    S(1+i)=N-I(1+i);
end
plot(t,I,t,S)
xlabel('时间')
ylabel('人数')
legend('患病者','易感者')
title('SI传染病模型')

亦可以简单理解为 $dIdt=βSI\frac{dI}{dt}=\beta S I$ ( $β\beta$ 是感染率)

SIS模型

此时只有易感者和患病者两类人群，但会反复发作的疾病。
患病者每天增加的比例：
$didt=i×(1−i)λ−μ⋅i\frac{di}{dt}=i\times(1-i)\lambda \color{red}{-\mu \cdot i}$
患病者每天增加的数量：
$dIdt=I×(1−I)λ/N−μ⋅I\frac{dI}{dt}=I \times (1-I)\lambda/N \color{red}{-\mu \cdot I}$
记 $α=λ/μ\alpha = \lambda / \mu$
模型可以解得：
当 $λ=μ\lambda = \mu$ 时， $i(t)=i0λti0+1i(t)=\frac{i_0}{\lambda t i_0 +1}$
当 $λ≠μ\lambda \not ={\mu}$ 时， $i(t)=11−α−1+i0(1−α−1)e(λ−μ)t(1−α−1)−i0i(t)=\frac{1}{1-\alpha^{-1}}+\frac{i_0(1-\alpha^{-1})e^{(\lambda - \mu)t}}{(1-\alpha^{-1})-i_0}$

clc;
close all;
clear all;
I=10;
N=10000;
S=N-I;
lemda=0.1;
mu=0.05;
t=1:365;
for i=1:(size(t,2)-1)
    I(1+i)=I(i)+I(i)*(N-I(i))*lemda/N-mu*I(i);
    S(1+i)=N-I(1+i);
end
plot(t,I,t,S)
xlabel('时间')
ylabel('人数')
legend('患病者','易感者')
title('SI传染病模型')

亦可以简单理解为 $dIdt=βIS−γI\frac{dI}{dt}=\beta I S-\gamma I$ ( $β是感染率γ是治愈率\beta 是感染率 \gamma 是治愈率$ )

SIR模型

患病后康复拥有永久性免疫，康复者可以理解为退出了传染链。
治愈率为 $μ\mu$ /天。


clc;
close all;
clear all;
I=10;
R=0;
N=10000;
S=N-I;
lemda=0.2;
mu=0.05;
t=1:365;
for i=1:(size(t,2)-1)
    I(1+i)=I(i)+I(i)*(N-I(i)-R(i))*lemda/N-mu*I(i);
    S(1+i)=S(i)-lemda*I(i)*S(i)/N;
    R(1+i)=N-I(1+i)-S(1+i);
end
plot(t,I,t,S,t,R)
xlabel('时间')
ylabel('人数')
legend('患病者','易感者','康复者')
title('SIR传染病模型')

存在微分方程：
$dSdt=−βIS\frac{dS}{dt}=- \beta I S$ ,
$dIdt=βIS−γI\frac{dI}{dt}=\beta I S - \gamma I$ ,
$dRdt=γI\frac{dR}{dt}=\gamma I$ ( $β是感染率γ是治愈率\beta 是感染率 \gamma 是治愈率$ )

SIRS SEIR

会反复感染反复患病，病毒存在潜伏期的复杂模型。

SIRS 康复者获得短暂的抗体
$dSdt=−βSI+αR\frac{dS}{dt}=-\beta S I+ \alpha R$ ,
$dIdt=betaSI−γI\frac{dI}{dt}=\\beta S I -\gamma I$ ,
$dRdt=γI−αR\frac{dR}{dt}=\gamma I -\alpha R$
( $β是感染率γ是治愈率α是复感率\beta 是感染率 \gamma 是治愈率 \alpha 是复感率$ )
SIER 先变成潜伏期患者再发病
$dSdt=−βIS\frac{dS}{dt}=- \beta I S$ ,
$dEdt=βIS−(α+γ1)E\frac{dE}{dt}=\beta I S -(\alpha +\gamma _1)E$
$dIdt=αE−γ2I\frac{dI}{dt}=\alpha E - \gamma _2 I$
$dRdt=γ1E+γ2I\frac{dR}{dt}=\gamma _1 E+\gamma _2 I$
$β\beta$ 是感染率， $γ1\gamma _1$ 为潜伏期康复率， $γ2\gamma _2$ 为患者康复率


function Message_Spread_Mode
tic
load 'Data\Link.txt';    %读入连接矩阵
% load '\Data\Point_X.txt'; %读入横坐标
% load '\Data\Point_Y.txt'; %读入纵坐标
%-------------------------------------------------------------------------%
%状态分布及状态转移概率SEIR
%0:易感状态S（Susceptible）  P_0_1; （P_0_3:预免疫系数）
%1:潜伏状态E（Exposed）      P_1_0；P_1_2；P_1_3
%2:染病状态I（Infected）     P_2_0；P_2_3
%3:免疫状态R（Recovered）    P_3_0
%-------------------------------------------------------------------------%
%计算各用户节点的度                                                          
De=sum(Link);                                                              %用户节点的度
%------------——————----参数设置与说明--------------------------------%
[M N]=size(Link);                                                          %连接矩阵的规模
I_E=0.6;                                                                   %潜伏期E用户的传染强度
I_I=0.9;                                                                   %发病期I用户的传染强度
lamda=sum(De)/M;                                                           %用户单位时间内平均发送信息的数量
%P_m1：用户预免疫系数
%State：用户所处状态State=zeros(1,M);0:表示易感状态（Susceptible）
%---------------------------------1---------------------------------------%
%先讨论用户预免疫系数P_m1对病毒传播的影响
TimeStep=50;%input('短信网络内病毒传播模拟时间：');
P_m1=[0.1,0.5,0.9];         %用户预免疫系数
% State=zeros(TimeStep,M);  %用户的状态  
G_t=5;                      %G_t：用户的免疫持续时间，反映了病毒的变异频率
F_t=5;                      %F_t：用户从发现病毒到杀毒并升级病毒库的时间
for i=1:length(P_m1)
    TimeLong_F=zeros(1,M); %用户处于染病期的时间长短
    TimeLong_E=zeros(1,M); %用户处于潜伏期的时间长短
    Sta=zeros(1,M);                                                      %用户的状态 
    %进行预免疫设定
    for j=1:M
        if rand(1)<=P_m1(i)
            Sta(j)=3;         %进入免疫状态
            TimeLong_E(j)=1;  %出入潜伏期的时间为1
        else
            continue;
        end
    end
    %状态转换
    %初始随机选择一个节点为病源点（此时不能选处于免疫状态的点）
    %问题：节点度大小存在差别，可能模拟出来的结果有出于
    %      为避免这个问题，我们取度最大的节点为病源节点，如果已免疫，则选次大的，一次下去
    [Number,Sta]=Select_Infected_Point(M,Sta,De);
    %Number：病源节点
    %State ：确定病源节点以后的节点状态矩阵
    State=zeros(TimeStep,M);
    Number_State=zeros(4,TimeStep);  %用户处于个状态的统计数量
    for t=1:TimeStep
        if t==1
            State(t,:)=Sta;
        else
            %模拟每个用户节点的状态
            for j=1:M
                %判断用户节点处于什么状态，然后根据其状态确定其转变情况
                if State(t-1,j)==0                          %此时处于易感状态0，可能向潜伏期转移
                    Num=Select_Number_Near(j,Link);         %找出节点j的邻居节点
                    P=zeros(1,length(Num));                 %邻居节点感染该节点的概率
                    for k=1:length(Num)
                        if State(t-1,Num(k))==1             %节点处于潜伏期E（1）
                            P(k)=I_E/De(Num(k))*sum((lamda.^([1:De(Num(k))]).*exp(-lamda))./...
                                (factorial([1:De(Num(k))]-1)));
                        else
                            if State(t-1,Num(k))==2          %节点处于染病期I（2）
                                P(k)=I_I/De(Num(k))*sum((lamda.^([1:De(Num(k))]).*exp(-lamda))./(factorial([1:De(Num(k))]-1)));
                            else
                                continue;
                            end
                        end
                    end
                    P_0_1=max(P);                       %节点感染病毒的概率
                    if rand<=P_0_1                      %此时节点进入潜伏期
                       State(t,j)=1;
                    else
                       State(t,j)=State(t-1,j); 
                    end
                else
                    if State(t-1,j)==1         %此时处于潜伏状态E，可能向易感S，染病I和免疫R转移
                        if rand<=1/(1+exp(-De(j)))                 %向染病状态I转移                
                            State(t,j)=2;
                            TimeLong_F(j)=TimeLong_F(j)+1;         %用户j处于染病状态的时间长短  
                        else
                            if rand<=1/(1+exp(-De(j)))             %向易感状态S转移           
                                State(t,j)=0;
                            else
                                if rand<=1/(1+exp(-De(j)))         %向免疫状态R转移
                                    State(t,j)=3;
                                    TimeLong_E(j)=TimeLong_E(j)+1; %免疫时间增加1
                                else
                                    State(t,j)=State(t-1,j);       %状态不变，依然为潜伏期E（1）
                                end
                            end
                        end
                    else
                        if State(t-1,j)==2        %此时处于欺染病状态I，可能向易感S，免疫R转移
                            if TimeLong_F(j)<=F_t         %表示此时用户不对病毒进行任何处理
                                State(t,j)=State(t-1,j);           %此时用户维持在原状态I
                                TimeLong_F(j)=TimeLong_F(j)+2;
                            else
                                %此时用户对进行杀毒并升级病毒库，进入免疫状态R
                                State(t,j)=3;
                                TimeLong_F(j)=0; %处于感染期（中毒状态）的时间长度
                                TimeLong_E(j)=1; %进入免疫期的时间长度
                            end
                        else
                            %此时用户处于免疫期
                            if TimeLong_E<=G_t   %病毒此时并未突变，维持原状态R（免疫状态）
                                State(t,j)=State(t-1,j);
                                TimeLong_E(j)=TimeLong_E(j)+1; %处于免疫期的时间增加
                            else
                                if rand<=1/G_t  %病毒突变，状态转移为易感状态S
                                    State(t,j)=0;
                                    TimeLong_E(j)=0;
                                else
                                    %此时用户状态依然不变
                                    State(t,j)=State(t-1,j);
                                    TimeLong_E(j)=TimeLong_E(j)+1; %处于免疫期的时间增加
                                end
                            end
                        end
                    end
                end
            end
        end
        %统计各状态的节点数量
        Number_State(1,t)=sum(State(t,:)==0);%处于易感状态S的总节点数量
        Number_State(2,t)=sum(State(t,:)==1);%处于易感状态E的总节点数量
        Number_State(3,t)=sum(State(t,:)==2);%处于易感状态I的总节点数量
        Number_State(4,t)=sum(State(t,:)==3);%处于易感状态R的总节点数量
        figure(i)
        if rem(t,3)==0
            plot([t-1 t],[Number_State(1,t-1) Number_State(1,t)],'md-'),hold on
            plot([t-1 t],[Number_State(2,t-1) Number_State(2,t)],'gh:'),hold on
            plot([t-1 t],[Number_State(3,t-1) Number_State(3,t)],'bs-.'),hold on
            plot([t-1 t],[Number_State(4,t-1) Number_State(4,t)],'k.-'),hold on
        else
            continue;
        end
        legend('易感状态Susceptible','潜伏状态Exposed','染病状态Infected','免疫状态Recovered')
        xlabel('模拟时间')
        ylabel('各状态的用户数量')
    end
end
P_m1=0.3;         %用户预免疫系数
% State=zeros(TimeStep,M);  %用户的状态  
% G_t=5;                      %G_t：用户的免疫持续时间，反映了病毒的变异频率
G_t=[1,5,9];
F_t=5;                      %F_t：用户从发现病毒到杀毒并升级病毒库的时间
for i=1:length(G_t)
    TimeLong_F=zeros(1,M); %用户处于染病期的时间长短
    TimeLong_E=zeros(1,M); %用户处于潜伏期的时间长短
    Sta=zeros(1,M);                                                      %用户的状态 
    %进行预免疫设定
    for j=1:M
        if rand(1)<=P_m1
            Sta(j)=3;         %进入免疫状态
            TimeLong_E(j)=1;  %出入潜伏期的时间为1
        else
            continue;
        end
    end
    %状态转换
    %初始随机选择一个节点为病源点（此时不能选处于免疫状态的点）
    %问题：节点度大小存在差别，可能模拟出来的结果有出于
    %      为避免这个问题，我们取度最大的节点为病源节点，如果已免疫，则选次大的，一次下去
    [Number,Sta]=Select_Infected_Point(M,Sta,De);
    %Number：病源节点
    %State ：确定病源节点以后的节点状态矩阵
    State=zeros(TimeStep,M);
    Number_State=zeros(4,TimeStep);  %用户处于个状态的统计数量
    for t=1:TimeStep
        if t==1
            State(t,:)=Sta;
        else
            %模拟每个用户节点的状态
            for j=1:M
                %判断用户节点处于什么状态，然后根据其状态确定其转变情况
                if State(t-1,j)==0                          %此时处于易感状态0，可能向潜伏期转移
                    Num=Select_Number_Near(j,Link);         %找出节点j的邻居节点
                    P=zeros(1,length(Num));                 %邻居节点感染该节点的概率
                    for k=1:length(Num)
                        if State(t-1,Num(k))==1             %节点处于潜伏期E（1）
                            P(k)=I_E/De(Num(k))*sum((lamda.^([1:De(Num(k))]).*exp(-lamda))./...
                                (factorial([1:De(Num(k))]-1)));
                        else
                            if State(t-1,Num(k))==2          %节点处于染病期I（2）
                                P(k)=I_I/De(Num(k))*sum((lamda.^([1:De(Num(k))]).*exp(-lamda))./...
                                    (factorial([1:De(Num(k))]-1)));
                            else
                                continue;
                            end
                        end
                    end
                    P_0_1=max(P);                       %节点感染病毒的概率
                    if rand<=P_0_1                      %此时节点进入潜伏期
                       State(t,j)=1;
                    else
                       State(t,j)=State(t-1,j); 
                    end
                else
                    if State(t-1,j)==1          %此时处于潜伏状态E，可能向易感S，染病I和免疫R转移
                        if rand<=1/(1+exp(-De(j)))                 %向染病状态I转移                
                            State(t,j)=2;
                            TimeLong_F(j)=TimeLong_F(j)+1;         %用户j处于染病状态的时间长短  
                        else
                            if rand<=1/(1+exp(-De(j)))             %向易感状态S转移           
                                State(t,j)=0;
                            else
                                if rand<=1/(1+exp(-De(j)))         %向免疫状态R转移
                                    State(t,j)=3;
                                    TimeLong_E(j)=TimeLong_E(j)+1; %免疫时间增加1
                                else
                                    State(t,j)=State(t-1,j);       %状态不变，依然为潜伏期E（1）
                                end
                            end
                        end
                    else
                        if State(t-1,j)==2           %此时处于欺染病状态I，可能向易感S，免疫R转移
                            if TimeLong_F(j)<=F_t          %表示此时用户不对病毒进行任何处理
                                State(t,j)=State(t-1,j);           %此时用户维持在原状态I
                                TimeLong_F(j)=TimeLong_F(j)+2;
                            else
                                %此时用户对进行杀毒并升级病毒库，进入免疫状态R
                                State(t,j)=3;
                                TimeLong_F(j)=0; %处于感染期（中毒状态）的时间长度
                                TimeLong_E(j)=1; %进入免疫期的时间长度
                            end
                        else
                            %此时用户处于免疫期
                            if TimeLong_E<=G_t(i)   %病毒此时并未突变，维持原状态R（免疫状态）
                                State(t,j)=State(t-1,j);
                                TimeLong_E(j)=TimeLong_E(j)+1; %处于免疫期的时间增加
                            else
                                if rand<=1/G_t(i)              %病毒突变，状态转移为易感状态S
                                    State(t,j)=0;
                                    TimeLong_E(j)=0;
                                else
                                    %此时用户状态依然不变
                                    State(t,j)=State(t-1,j);
                                    TimeLong_E(j)=TimeLong_E(j)+1; %处于免疫期的时间增加
                                end
                            end
                        end
                    end
                end
            end
        end
        %统计各状态的节点数量
        Number_State(1,t)=sum(State(t,:)==0);%处于易感状态S的总节点数量
        Number_State(2,t)=sum(State(t,:)==1);%处于易感状态E的总节点数量
        Number_State(3,t)=sum(State(t,:)==2);%处于易感状态I的总节点数量
        Number_State(4,t)=sum(State(t,:)==3);%处于易感状态R的总节点数量
        figure(i+5)
        if rem(t,3)==0
            plot([t-1 t],[Number_State(1,t-1) Number_State(1,t)],'md-'),hold on
            plot([t-1 t],[Number_State(2,t-1) Number_State(2,t)],'gh:'),hold on
            plot([t-1 t],[Number_State(3,t-1) Number_State(3,t)],'bs-.'),hold on
            plot([t-1 t],[Number_State(4,t-1) Number_State(4,t)],'k.-'),hold on
        else
            continue;
        end
        legend('易感状态Susceptible','潜伏状态Exposed','染病状态Infected','免疫状态Recovered')
        xlabel('模拟时间')
        ylabel('各状态的用户数量')
    end
end
P_m1=0.3;                   %用户预免疫系数
% State=zeros(TimeStep,M);  %用户的状态  
% G_t=5;                      %G_t：用户的免疫持续时间，反映了病毒的变异频率
G_t=5;
F_t=[1,5,9];                        %F_t：用户从发现病毒到杀毒并升级病毒库的时间
for i=1:length(F_t)
    TimeLong_F=zeros(1,M); %用户处于染病期的时间长短
    TimeLong_E=zeros(1,M); %用户处于潜伏期的时间长短
    Sta=zeros(1,M);                                                      %用户的状态 
    %进行预免疫设定
    for j=1:M
        if rand(1)<=P_m1
            Sta(j)=3;         %进入免疫状态
            TimeLong_E(j)=1;  %出入潜伏期的时间为1
        else
            continue;
        end
    end
    %状态转换
    %初始随机选择一个节点为病源点（此时不能选处于免疫状态的点）
    %问题：节点度大小存在差别，可能模拟出来的结果有出于
    %      为避免这个问题，我们取度最大的节点为病源节点，如果已免疫，则选次大的，一次下去
    [Number,Sta]=Select_Infected_Point(M,Sta,De);
    %Number：病源节点
    %State ：确定病源节点以后的节点状态矩阵
    State=zeros(TimeStep,M);
    Number_State=zeros(4,TimeStep);  %用户处于个状态的统计数量
    for t=1:TimeStep
        if t==1
            State(t,:)=Sta;
        else
            %模拟每个用户节点的状态
            for j=1:M
                %判断用户节点处于什么状态，然后根据其状态确定其转变情况
                if State(t-1,j)==0                          %此时处于易感状态0，可能向潜伏期转移
                    Num=Select_Number_Near(j,Link);         %找出节点j的邻居节点
                    P=zeros(1,length(Num));                 %邻居节点感染该节点的概率
                    for k=1:length(Num)
                        if State(t-1,Num(k))==1             %节点处于潜伏期E（1）
                            P(k)=I_E/De(Num(k))*sum((lamda.^([1:De(Num(k))]).*exp(-lamda))./...
                                (factorial([1:De(Num(k))]-1)));
                        else
                            if State(t-1,Num(k))==2          %节点处于染病期I（2）
                                P(k)=I_I/De(Num(k))*sum((lamda.^([1:De(Num(k))]).*exp(-lamda))./...
                                    (factorial([1:De(Num(k))]-1)));
                            else
                                continue;
                            end
                        end
                    end
                    P_0_1=max(P);                       %节点感染病毒的概率
                    if rand<=P_0_1                      %此时节点进入潜伏期
                       State(t,j)=1;
                    else
                       State(t,j)=State(t-1,j); 
                    end
                else
                    if State(t-1,j)==1      %此时处于潜伏状态E，可能向易感S，染病I和免疫R转移
                        if rand<=1/(1+exp(-De(j)))                 %向染病状态I转移                
                            State(t,j)=2;
                            TimeLong_F(j)=TimeLong_F(j)+1;         %用户j处于染病状态的时间长短  
                        else
                            if rand<=1/(1+exp(-De(j)))             %向易感状态S转移           
                                State(t,j)=0;
                            else
                                if rand<=1/(1+exp(-De(j)))         %向免疫状态R转移
                                    State(t,j)=3;
                                    TimeLong_E(j)=TimeLong_E(j)+1; %免疫时间增加1
                                else
                                    State(t,j)=State(t-1,j);       %状态不变，依然为潜伏期E（1）
                                end
                            end
                        end
                    else
                        if State(t-1,j)==2           %此时处于欺染病状态I，可能向易感S，免疫R转移
                            if TimeLong_F(j)<=F_t(i)   %表示此时用户不对病毒进行任何处理
                                State(t,j)=State(t-1,j);           %此时用户维持在原状态I
                                TimeLong_F(j)=TimeLong_F(j)+2;
                            else
                                %此时用户对进行杀毒并升级病毒库，进入免疫状态R
                                State(t,j)=3;
                                TimeLong_F(j)=0; %处于感染期（中毒状态）的时间长度
                                TimeLong_E(j)=1; %进入免疫期的时间长度
                            end
                        else
                            %此时用户处于免疫期
                            if TimeLong_E<=G_t          %病毒此时并未突变，维持原状态R（免疫状态）
                                State(t,j)=State(t-1,j);
                                TimeLong_E(j)=TimeLong_E(j)+1; %处于免疫期的时间增加
                            else
                                if rand<=1/G_t              %病毒突变，状态转移为易感状态S
                                    State(t,j)=0;
                                    TimeLong_E(j)=0;
                                else
                                    %此时用户状态依然不变
                                    State(t,j)=State(t-1,j);
                                    TimeLong_E(j)=TimeLong_E(j)+1; %处于免疫期的时间增加
                                end
                            end
                        end
                    end
                end
            end
        end
        %统计各状态的节点数量
        Number_State(1,t)=sum(State(t,:)==0);%处于易感状态S的总节点数量
        Number_State(2,t)=sum(State(t,:)==1);%处于易感状态E的总节点数量
        Number_State(3,t)=sum(State(t,:)==2);%处于易感状态I的总节点数量
        Number_State(4,t)=sum(State(t,:)==3);%处于易感状态R的总节点数量
        figure(i+10)
        if rem(t,3)==0
            plot([t-1 t],[Number_State(1,t-1) Number_State(1,t)],'md-'),hold on
            plot([t-1 t],[Number_State(2,t-1) Number_State(2,t)],'gh:'),hold on
            plot([t-1 t],[Number_State(3,t-1) Number_State(3,t)],'bs-.'),hold on
            plot([t-1 t],[Number_State(4,t-1) Number_State(4,t)],'k.-'),hold on
        else
            continue;
        end
        legend('易感状态Susceptible','潜伏状态Exposed','染病状态Infected','免疫状态Recovered')
        xlabel('模拟时间')
        ylabel('各状态的人口数量')
    end
end
toc