疫情传播SEIR模型(python)

假设有N个人,两两之间有接触的关系为p。

在这里插入图片描述

事件开始: 选中一个人变为潜伏者 E E E

时间推移:

  • 健康者 S S S:若与 x x x个感染者 I I I接触,则有 1 − ( 1 − w ) x 1-(1-w)^x 1(1w)x的概率变为潜伏者 E E E
  • 潜伏者 E E E:有 β \beta β的概率病发变为感染者 I I I,有 μ \mu μ的概率变为终态 R R R
  • 感染者 I I I:有 μ \mu μ的概率变为终态 R R R

参数设置:

N = 1000 , p = 0.006 w = 0.2 , β = 0.5 , μ = 0.2 N = 1000,p = 0.006\\ w=0.2,\beta=0.5,\mu=0.2 N=1000,p=0.006w=0.2,β=0.5,μ=0.2

代码:

import numpy as np
import matplotlib.pyplot as plt

def random_rate(p):
    return np.random.rand() < p


class SEIR():
    def __init__(self):
        self.N = int(1000)
        self.p = 0.006
        self.w = 0.2
        self.b = 0.5
        self.u = 0.2
        self.t = 100
        self.link = []
        for i in range(self.N):
            self.link.append([])
        # 0 S  1 E  2 I  3 R
        self.category = np.array(np.zeros([self.N], dtype=int))
        self.category[np.random.randint(self.N)] = 2

    def generate_graph(self):
        for i in range(self.N):
            for j in range(i + 1, self.N):
                if random_rate(self.p):
                    self.link[i].append(j)
                    self.link[j].append(i)

    def print_graph(self):
        # for list in self.link:
        #     print(list)

        mxdeg=0
        for i in range(self.N):
            mxdeg=max(mxdeg,self.link[i].__len__())
        mxdeg+=1

        deg=[0]*mxdeg
        for i in range(self.N):
            deg[self.link[i].__len__()]+=1
        for i in range(mxdeg):
            deg[i]/=1.0*self.N
        plt.plot(np.linspace(0,mxdeg-1,mxdeg),deg)
        plt.xlabel('degree')
        plt.ylabel('P')
        plt.title('The probability distribution of degree')
        plt.show()

    def time_advancing(self):
        probability = np.array(np.zeros([self.N], dtype=float))
        for i in range(self.N):
            if self.category[i] != 0:
                continue
            count = 0
            for u in self.link[i]:
                if self.category[u] == 2:
                    count += 1
            probability[i] = 1 - np.float_power(1 - self.w, count)
        for i in range(self.N):
            if self.category[i] == 0:
                if random_rate(probability[i]):
                    self.category[i] = 1
            elif self.category[i] == 1:
                p = np.random.rand()
                if p < self.b:
                    self.category[i] = 2
                elif p < self.b + self.u:
                    self.category[i] = 3
            elif self.category[i] == 2:
                if random_rate(self.u):
                    self.category[i] = 3

    def count_number(self):
        res=[0,0,0,0]
        for i in range(self.N):
            res[self.category[i]]+=1
        return res

    def running(self):
        result=np.matrix(np.zeros([4,self.t+1], dtype=int))
        for i in range(self.t+1):
            numbers=self.count_number()
            for j in range(4):
                result[j,i]=numbers[j]
            self.time_advancing()
        self.show_result(result)

    def show_result(self,result):
        plt.title('SEIR model result')
        plt.xlabel('time')
        plt.ylabel('person')
        x=np.linspace(0,self.t,self.t+1)
        s1='SEIR'
        s2='>v<^'
        for idx in range(4):
            y=[]
            for i in range(self.t+1):
                y.append(result[idx,i])
            plt.plot(x,y,'-{}'.format(s2[idx]),label=s1[idx])
        plt.legend()
        plt.show()

model = SEIR()
model.generate_graph()
model.print_graph()
model.running()

初始图模型:
在这里插入图片描述
时间推移:
在这里插入图片描述

已标记关键词 清除标记
相关推荐
©️2020 CSDN 皮肤主题: 1024 设计师:白松林 返回首页