MENU

Python实现传染病动力模型计算

• February 13, 2020 • Read: 5823 • Knowledge,Python,Numerical Analysis

$SI$模型

$N$表示一个城市的总人口,在不考虑出生人口和死亡人口的情况下,$I$代表已经感染人的个数,称为患者数(或者僵尸数);这一模型将人群分为两类,一类是感染的人,另外一类是未被感染的人,也称为易感染人,数量用$S$表示。自然可以有$N = I + S$。

每一天一个患者能够传染的速率为$\lambda$,那么已有的感染人的数量乘以这个感染传播速率就得到新增的感染人数,$I \times \lambda$。这些都是对于每一天来说的,需要注意的是 前一天的情况会得到下一天的情况,因此,这一个过程是一种完全的迭代实现。

令感染者在人群中的比例为$i = \frac{I}{N}$,易感染者在人群中的比例为$s = \frac{S}{N}$,因为对于人群的划分只分为两类,所以得到$1 = i + s$。

依据上述的表达,可以计算出每一天新增的感染人数为$\lambda \times iN \times s$,这里乘以$s$是因为,只要感染者增加,那么易感人数必定减少,而新增的感染者一定是从未感染者中来的。这样新增感染者比例的增加量为$\lambda i s$。同时借助于$1 = s + i$,得到$s = 1- i$,所以后一天感染者所占的比例为$i + \lambda i (1 - i)$。由此通过该迭代关系可以完成预测。

以一个例子实现,假设总人口为三百七万人(3.7e6),初始时刻感染的人为10个,感染速率$\lambda$设定为1,即每一天一个感染者会感染一个人,模型总时间为100天。这里不讨论参数如何确定,只是单纯地为了计算,结果并不能代表现实中的情况。

import numpy as np
import matplotlib.pyplot as plt

class SI(object):
    
    def __init__(self, N, T, m, alpha):
        self.N = N
        self.T = T
        self.m = m
        self.alpha = alpha
        self.get_infective_ratio()
        
    def get_infective_ratio(self):
        ii = np.zeros(self.T)
        ii[0] = self.m / self.N
        for t in np.arange(self.T - 1):
            ii[t+1] = ii[t] + self.alpha *  ii[t] * (1 - ii[t])
        self.s = 1 - ii
        self.i = ii
        return self.i, self.s
    
    def visulization(self):
        plt.style.use('ggplot')
        fig, ax = plt.subplots()
        ax.plot(self.i, color='tab:red', linewidth=2, label='I')
        ax.plot(self.s, color='tab:blue', linewidth=2, label='S')
        #ax.grid(True)
        ax.set_xlabel('Days/ day')
        ax.set_ylabel('Ratio')
        ax.set_title('SI Model')
        ax.legend(shadow=True,loc='best')
        fig.savefig('si.png', dpi=600)
        plt.show()

demo = SI(3.7e6, 100, 10, 1)
demo.visulization()

si-model

$SI$模型是一个最简单的模型,这一模型中人群分为感染者和未感染者两种,因为总人口不变,所以这两者之间的关系一定是“此消彼长”的,感染者数量的增加一定会导致未感染者数量的下降,所以,直观地可以意识到,感染者所占的比例已开始会缓慢增加,随后会随着已有感染者数量的速度加快,但是这不是无限的,因为未感染者数量在持续的下降,直到两者数量相等时,出现分水岭,感染比例速度开始下降。这一模型过于理想化,考虑的因素也过于简单。

$SIR$模型

$SI$模型中没有考虑死亡和治愈的人的影响,实际情况中,死亡的人不会再进行疾病的传播,如果患者被治愈其体内会有抗体存在,应当来说是不会再被感染的,这两种情况造成的影响都是相同的,都是“去除”的效果,这也是$R$的意义(remove)。

$SIR$模型是在$SI$模型上的一个改进,基本的原理不变,只是加入了$R$的影响。移出者$R$一定是来自于已有的感染人群,不论是死亡还是治愈,设定移出率为$\gamma$。这一模型中人群被分为了3类。

$S$易感染者:由于疾病一直在传播,所以易感染的人数量一直在下降,它的变化数量就是增加的感染者数量,$\lambda \times N \times i \times s$,则对应的变化比例为$\lambda i s$。

$I$感染者:不考虑移除为$\lambda N i s$,现在加入了移除$R$,那么需要减去$R$的数量,得到$\lambda N i s - \gamma N i$,得到对应的比例为$\lambda is - \gamma i$。

$R$移出者:移出者是来自于已经感染的人群,即$\gamma N i$,所占的比例为$\gamma i$。

这样基于上述的描述,利用迭代计算可以计算出设定的天数内各个人群比例的变化值。

为了演示计算,设定总人口依然为三百七十万人,初始时刻的感染者为6人,$\lambda=1$,Remove速率为0.3,分析的总时间是100天。

class SIR(object):
    
    def __init__(self, N, T, m, alpha, gamma):
        self.N = N
        self.T = T
        self.m = m
        self.alpha = alpha
        self.gamma = gamma
        self.get_result()
    
    def get_result(self):
        i0 = np.zeros(self.T)
        s0 = np.zeros(self.T)
        r0 = np.zeros(self.T)
        i0[0] = self.m / self.N
        s0[0] = 1 - i0[0]
        for t in np.arange(self.T-1):
            i0[t+1] = i0[t] + self.alpha * i0[t] * s0[t] - i0[t] * self.gamma
            s0[t+1] = s0[t] - self.alpha * i0[t] * s0[t]
            r0[t+1] = r0[t] + self.gamma * i0[t]
        self.i = i0
        self.s = s0
        self.r = r0
        return self.i, self.s, self.r
    
    def visulization(self):
        plt.style.use('ggplot')
        fig, ax = plt.subplots()
        ax.plot(self.i, color='tab:red', linewidth=2, label='I')
        ax.plot(self.s, color='tab:blue', linewidth=2, label='S')
        ax.plot(self.r, color='tab:orange', linewidth=2, label='R')
        #ax.grid(True)
        ax.set_xlabel('Days/ day')
        ax.set_ylabel('Ratio')
        ax.set_title('SIR Model')
        ax.legend(shadow=True,loc='best')
        fig.savefig('SIR.png', dpi=600)
        plt.show()
        
demo1 = SIR(3.7e6, 100, 6, 1, 0.3)
demo1.visulization()

sir-model

这里只是理想地进行一个计算,具体情况影响因素很多,上述的模型都是过于理想化的,但是能反映出感染率等因素对于疾病传播规模的影响,其中$\lambda$值在实际情况中就体现为防控措施。

Last Modified: February 16, 2020
Archives Tip
QR Code for this page
Tipping QR Code