Python小白的数学建模课-B3. 新冠疫情 SIS模型

Posted 2021-06-26 youcans

传染病的数学模型是数学建模中的典型问题，常见的传染病模型有 SI、SIR、SIRS、SEIR 模型。

SIS 模型型将人群分为 S 类和 I 类，考虑患病者可以治愈而变成易感者，但不考虑免疫期。

本文详细给出了 SIS 模型的建模、例程、运行结果和模型分析，让小白都能懂。

『Python小白的数学建模课 @ Youcans』带你从数模小白成为国赛达人。

---

1. 疫情传播 SIS 模型

传染病动力学是对传染病进行定量研究的重要方法。它依据种群繁衍迁移的特性、传染病在种群内产生及传播的机制、医疗与防控条件等外部因素，建立可以描述传染病动力学行为的数学模型，通过对模型进行定性、定量分析和数值计算，模拟传染病的传播过程，预测传染病的发展趋势，研究防控策略的作用。

1.1 SI 模型

SI 模型把人群分为易感者（S类）和患病者（I类）两类，易感者（S类）与患病者（I类）有效接触即被感染，变为患病者，无潜伏期、无治愈情况、无免疫力。

SI 模型适用于只有易感者和患病者两类人群，且无法治愈的疾病。

按照 SI 模型，最终所有人都会被传染而变成病人，这是因为模型中没有考虑病人可以治愈。因此只能是健康人患病，而患病者不能恢复健康（甚至也不会死亡，而是不断传播疫情），所以终将全部被传染。

欢迎关注 『Python小白的数学建模课 @ Youcans』 系列，持续更新

1.2 SIS 模型

SIS 模型将人群分为 S 类和 I 类，考虑患病者（I 类）可以治愈而变成易感者（S 类），但不考虑免疫期，因此患病者（I 类）治愈变成易感者以后还可以被感染而变成患病者。

SIS 模型适用于只有易感者和患病者两类人群，可以治愈，但会反复发作的疾病，例如脑炎、细菌性痢疾等治愈后也不具有免疫力的传染病。

SIS 模型假设：

1. 考察地区的总人数 N 不变，即不考虑生死或人口流动；
2. 人群分为易感者（S类）和患病者（I类）两类；
3. 易感者（S类）与患病者（I类）有效接触即被感染，变为患病者；患病者（I类）可被治愈而变为易感者，无潜伏期、无免疫力；
4. 每个患病者每天有效接触的易感者的平均人数（日接触数）是 $\lambda$，称为日接触率；
5. 每天被治愈的患病者人数占患病者总数的比例为 $\mu$ ，即日治愈率；
6. 将第 t 天时 S类、I 类人群的占比记为 $s(t)$、$i(t)$，数量为 $S(t)$、$I(t)$；初始日期 $t=0$ 时， S类、I 类人群占比的初值为 $s_0$、$i_0$。

需要说明的是，不考虑生死或人口流动，通常是由于考虑一个封闭环境而且假定疫情随时间的变化比生死、迁移随时间的变化显著得多, 因此后者可以忽略不计。

SIS 模型的微分方程：

由
$$N\frac{di}{dt}=N\lambda si-N\mu i$$
得：
$$\frac{di}{dt}=\lambda i(1-i)-\mu i，i(0)=i_0$$
由日治愈率 $\mu$ 可知平均治愈天数为 $1/\mu$，也称平均传染期。定义 $\sigma =\lambda /\mu$，其含义是每个病人在传染期内所传染的平均人数，称为传染期接触数。例如，平均传染期 $1/\mu =5$，日接触率 $\lambda =2$（每天传染 2人），则传染期接触数 $\sigma =10$。

SIS 模型的解析解为：
$$\{\begin{array}{l}\begin{array}{rlr}& i(t)=\frac{i_0}{1+\lambda t{i}_0}& ,\lambda =\mu \\ & i(t)=[\frac{\lambda }{\lambda -\mu }+(\frac{1}{i_0}-\frac{\lambda }{\lambda -\mu })\ast e^{-(\lambda -\mu )t}{]}^{-1}& ,\lambda \ne \mu \end{array}\end{array}$$

注意：网上有些博文中解析解的公式误写成 $exp((\lambda -\mu )t)$ ，漏掉了一个负号。

---

2. SIS 模型的 Python 编程

2.1 Scipy 工具包求解 SIS 模型

SIS 模型是常微分方程初值问题，可以使用 Scipy 工具包的 scipy.integrate.odeint() 函数求数值解。

scipy.integrate.odeint(func, y0, t, args=())

**scipy.integrate.odeint() **是求解微分方程的具体方法，通过数值积分来求解常微分方程组。

odeint() 的主要参数：

• func: callable(y, t, …) 　　导数函数 $f(y,t)$ ，即 y 在 t 处的导数，以函数的形式表示
• y0: array：　　初始条件 $y_0$，对于常微分方程组 $y_0$ 则为数组向量
• t: array：　　求解函数值对应的时间点的序列。序列的第一个元素是与初始条件 $y_0$ 对应的初始时间 $t_0$；时间序列必须是单调递增或单调递减的，允许重复值。
• args: 向导数函数 func 传递参数。当导数函数 $f(y,t,p1,p2,..)$ 包括可变参数 p1,p2… 时，通过 args =(p1,p2,…) 可以将参数p1,p2… 传递给导数函数 func。

odeint() 的返回值：

• y: array 　　数组，形状为 (len(t),len(y0)，给出时间序列 t 中每个时刻的 y 值。

odeint() 的编程步骤：

1. 导入 scipy、numpy、matplotlib 包；
2. 定义导数函数 $f(i,t)=\lambda i(1-i)-\mu i$ ；
3. 定义初值 $y_0$ 和 $y$ 的定义区间 $[t_0,t]$；
4. 调用 odeint() 求 $y$ 在定义区间 $[t_0,t]$ 的数值解。

2.2 Python例程：SIS 模型的解析解与数值解

# 1. SIS 模型，常微分方程，解析解与数值解的比较
from scipy.integrate import odeint  # 导入 scipy.integrate 模块
import numpy as np  # 导入 numpy包
import matplotlib.pyplot as plt  # 导入 matplotlib包

def dy_dt(y, t, lamda, mu):  # SIS 模型，导数函数
    dy_dt = lamda*y*(1-y) - mu*y  # di/dt = lamda*i*(1-i)-mu*i
    return dy_dt

# 设置模型参数
number = 1e5  # 总人数
lamda = 1.2  # 日接触率, 患病者每天有效接触的易感者的平均人数
sigma = 2.5  # 传染期接触数
mu = lamda/sigma  # 日治愈率, 每天被治愈的患病者人数占患病者总数的比例
fsig = 1-1/sigma
y0 = i0 = 1e-5  # 患病者比例的初值
tEnd = 50  # 预测日期长度
t = np.arange(0.0,tEnd,1)  # (start,stop,step)
print("lamda={}\\tmu={}\\tsigma={}\\t(1-1/sig)={}".format(lamda,mu,sigma,fsig))

# 解析解
if lamda == mu:
    yAnaly = 1.0/(lamda*t +1.0/i0)
else:
    yAnaly= 1.0/((lamda/(lamda-mu)) + ((1/i0)-(lamda/(lamda-mu))) * np.exp(-(lamda-mu)*t))
# odeint 数值解，求解微分方程初值问题
ySI = odeint(dy_dt, y0, t, args=(lamda,0))  # SI 模型
ySIS = odeint(dy_dt, y0, t, args=(lamda,mu))  # SIS 模型

# 绘图
plt.plot(t, yAnaly, '-ob', label='analytic')
plt.plot(t, ySIS, ':.r', label='ySIS')
plt.plot(t, ySI, '-g', label='ySI')

plt.title("Comparison between analytic and numerical solutions")
plt.axhline(y=fsig,ls="--",c='c')  # 添加水平直线
plt.legend(loc='best')  # youcans
plt.axis以上是关于Python小白的数学建模课-B3. 新冠疫情 SIS模型的主要内容，如果未能解决你的问题，请参考以下文章
 
 Python小白的数学建模课-A3. 12个新冠疫情数模竞赛赛题与点评
 Python小白的数学建模课-A3. 10个新冠疫情数模竞赛赛题与点评
 Python小白的数学建模课-B2. 新冠疫情 SI模型
 Python小白的数学建模课-B5. 新冠疫情 SEIR模型
 Python小白的数学建模课-B5. 新冠疫情 SEIR模型
 Python小白的数学建模课-B6. 新冠疫情 SEIR 改进模型