Python小白的數學建模課-B3. 新冠疫情 SIS模型

時間 2021-07-23

標籤 html python 算法編程數組函數工具 spa code orm 欄目 Python 简体版

原文原文鏈接

傳染病的數學模型是數學建模中的典型問題，常見的傳染病模型有 SI、SIR、SIRS、SEIR 模型。html

SIS 模型型將人羣分爲 S 類和 I 類，考慮患病者能夠治癒而變成易感者，但不考慮免疫期。python

本文詳細給出了 SIS 模型的建模、例程、運行結果和模型分析，讓小白都能懂。算法

『Python小白的數學建模課 @ Youcans』帶你從數模小白成爲國賽達人。編程

Python小白的數學建模課-09.微分方程模型
 Python小白的數學建模課-B2.新冠疫情 SI模型
 Python小白的數學建模課-B3.新冠疫情 SIS模型
 Python小白的數學建模課-B4.新冠疫情 SIR模型
 Python小白的數學建模課-B5.新冠疫情 SEIR模型
 Python小白的數學建模課-B6.改進 SEIR疫情模型數組

1. 疫情傳播 SIS 模型

傳染病動力學是對傳染病進行定量研究的重要方法。它依據種羣繁衍遷移的特性、傳染病在種羣內產生及傳播的機制、醫療與防控條件等外部因素，創建能夠描述傳染病動力學行爲的數學模型，經過對模型進行定性、定量分析和數值計算，模擬傳染病的傳播過程，預測傳染病的發展趨勢，研究防控策略的做用。函數

1.1 SI 模型

SI 模型把人羣分爲易感者（S類）和患病者（I類）兩類，易感者（S類）與患病者（I類）有效接觸即被感染，變爲患病者，無潛伏期、無治癒狀況、無免疫力。工具

SI 模型適用於只有易感者和患病者兩類人羣，且沒法治癒的疾病。spa

按照 SI 模型，最終全部人都會被傳染而變成病人，這是由於模型中沒有考慮病人能夠治癒。所以只能是健康人患病，而患病者不能恢復健康（甚至也不會死亡，而是不斷傳播疫情），因此終將所有被傳染。code

1.2 SIS 模型

SIS 模型將人羣分爲 S 類和 I 類，考慮患病者（I 類）能夠治癒而變成易感者（S 類），但不考慮免疫期，所以患病者（I 類）治癒變成易感者之後還能夠被感染而變成患病者。orm

SIS 模型適用於只有易感者和患病者兩類人羣，能夠治癒，但會反覆發做的疾病，例如腦炎、細菌性痢疾等治癒後也不具備免疫力的傳染病。

SIS 模型假設：

考察地區的總人數 N 不變，即不考慮生死或人口流動；
人羣分爲易感者（S類）和患病者（I類）兩類；
易感者（S類）與患病者（I類）有效接觸即被感染，變爲患病者；患病者（I類）可被治癒而變爲易感者，無潛伏期、無免疫力；
每一個患病者天天有效接觸的易感者的平均人數（日接觸數）是 $\lambda$，稱爲日接觸率；
天天被治癒的患病者人數佔患病者總數的比例爲 $\mu$ ，即日治癒率；
將第 t 天時 S類、I 類人羣的佔比記爲 $s(t)$、$i(t)$，數量爲 $S(t)$、$I(t)$；初始日期 $t=0$ 時， S類、I 類人羣佔比的初值爲 $s_0$、$i_0$。

須要說明的是，不考慮生死或人口流動，一般是因爲考慮一個封閉環境並且假定疫情隨時間的變化比生死、遷移隨時間的變化顯著得多, 所以後者能夠忽略不計。

SIS 模型的微分方程：

由

\[\begin{align*} N\frac{di}{dt} = N \lambda s i - N \mu i \end{align*} \]

得：

\[\begin{align*} \frac{di}{dt} = \lambda i (1-i) - \mu i，\ i(0) = i_0 \end{align*} \]

由日治癒率 $\mu$ 可知平均治癒天數爲 $1/\mu$，也稱平均傳染期。定義 $\sigma = \lambda / \mu$，其含義是每一個病人在傳染期內所傳染的平均人數，稱爲傳染期接觸數。例如，平均傳染期 $1/\mu = 5$，日接觸率 $\lambda = 2$（天天傳染 2人），則傳染期接觸數 $\sigma = 10$。

SIS 模型的解析解爲：

\[\begin{cases} \begin{aligned} & 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})*e^{-(\lambda - \mu) t}]^{-1} &,\lambda \neq \mu\\ \end{aligned} \end{cases}\\ \]

注意：網上有些博文中解析解的公式誤寫成 $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() 的編程步驟：

導入 scipy、numpy、matplotlib 包；
定義導數函數 $f(i,t)=\lambda i (1-i)- \mu i$ ；
定義初值 $y_0$ 和 $y$ 的定義區間 $[t_0,\ t]$；
調用 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([0, 50, -0.1, 1.1])
plt.show()

2.3 SIS 模型解析解與數值解的比較

本圖爲例程 2.2 的運行結果，圖中對解析解（藍色）與使用 odeint() 獲得的數值解（紅色）進行比較。在該例中，沒法觀察到解析解與數值解的差別，代表數值解的偏差很小。

本圖也比較了對相同日接觸率和患病者初值下 SI模型與 SIS模型進行了比較。SI 模型更早進入爆發期，最終收斂到 100%；SIS 模型下進入爆發期較晚，患病者的比例最終收斂到某個常數（與模型參數有關）。

考察 SI 模型與 SIS模型的關係，顯然 SI 模型是 SIS 模型在 $\mu = 0$ 時的特殊狀況。

3. SIS 模型參數的影響

對於 SIS 模型，須要考慮日接觸率 $\lambda$ 與日治癒率 $\mu$ 的關係、患病者比例的初值 $i_0$ 的影響，總人數 N 沒有影響。

3.1 日接觸率 $\lambda$ 與日治癒率 $\mu$ 關係的影響

直觀地考慮，若是天天治癒的人數高於感染的人數，則疫情逐漸好轉，不然疫情逐漸嚴重。所以日接觸率 $\lambda$ 與日治癒率 $\mu$ 的關係很是關鍵，這就是傳染期接觸數 $\sigma = \lambda / \mu$ 的意義。

（1） $\sigma \leq 1$

當 $\sigma<1$ 時，傳染期接觸數小於 1，日接觸率小於日治癒率，患病率單調降低，最終清零，與患病率初值無關。 $\sigma$ 越小，疫情清零速度越快； $\sigma$ 越接近於 1，疫情清零越慢，但最終仍將清零。

分析其實際意義，傳染期接觸數小於 1，代表在傳染期內通過接觸而使易感者變成患病者的數量，小於在傳染期內治癒的患病者的數量，所以患病者數量、比例都會逐漸下降，因此最終能夠清零，稱爲無病平衡點。

當 $\sigma=1$ 時，不論患病率初值如何，患病率也是單調降低，最終趨近於 0。雖然在數學上患病率只能趨近於 0 而不等於 0，但考慮到總人數 N 是有限的，而患病者和易感者人數須要取整，所以 $\sigma=1$ 時最終也會清零。

（2） $\sigma > 1$

當 $\sigma>1$ 時，傳染期接觸數大於 1，日接觸率大於日治癒率，患病率的升降有兩種狀況：

當患病率很低時，患病者人數少而易感者人數多，患病率上升；但隨着患病率增大，患病者愈來愈多而易感者愈來愈少，患病率雖然仍然上升但上升速度趨緩，最終趨於定值。
當患病率很高時，患病者人數多而易感者人數少，患病率降低；但隨着患病率減少，患病者愈來愈少而易感者愈來愈多，患病率雖然仍然降低但降低速度趨緩，最終也趨於相同的定值。
患病率最終都會收斂到穩態特徵值 $i_\infty=1-1/\sigma$。當 $i_0>i_\infty$ 即患病率初值大於穩態特徵值時，疫情曲線單調上升收斂；當 $i_0<i_\infty$ 即患病率初值小於穩態特徵值時，疫情曲線單調降低收斂；當 $i_0 = i_\infty$ 時，患病率始終大於穩態特徵值，疫情曲線爲水平直線。

這代表，當 $\sigma>1$ 時疫情終將穩定但不會清零，而是長期保持必定的患病率，稱爲地方病平衡點。

當 $\sigma=1$ 時，不論患病率初值如何，患病率都單調降低並最終趨於 0。

3.2 傳染期接觸數 $\sigma$ 與 $ di/dt$ 的關係

患病率的一階導數 $di/dt$ 的變化曲線，代表不論傳染期接觸數和初值如何，患病率的變化率都將收斂到 0，所以疫情終將穩定。當 $\sigma<1$ 時， $di/dt$ 始終是負值，單調上升趨近於 0；當 $\sigma>1$ 時， $di/dt$ 始終是正值，先上升達到峯值後再逐漸減少趨近於 0。

本圖爲患病率 $i(t)$ 與一階導數 $di/dt$ 在不一樣傳染期接觸數下的關係曲線（相空間圖）。當 $\sigma\leq 1$ 時，曲線收斂到原點 $(0,0)$，即存在無病平衡點；當 $\sigma>1$ 時，曲線收斂到 $(1-1/\sigma,0)$，即存在地方病平衡點。

3.3 Python例程：傳染期接觸數 $\sigma$ 與 $ di/dt$ 的關係

# 4. SIS 模型，模型參數對 di/dt的影響
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 = np.array((0.1, 0.5, 0.8, 0.95, 1.0))  # 傳染期接觸數
sigma = np.array((0.5, 0.8, 1.0, 1.5, 2.0, 3.0))  # 傳染期接觸數
y0 = i0 = 0.05  # 患病者比例的初值
tEnd = 100  # 預測日期長度
t = np.arange(0.0,tEnd,0.1)  # (start,stop,step)

for p in sigma:
    ySIS = odeint(dy_dt, y0, t, args=(lamda,lamda/p))  # SIS 模型
    yDeriv = lamda*ySIS*(1-ySIS) - ySIS*lamda/p
    # plt.plot(t, yDeriv, '-', label=r"$\sigma$ = {}".format(p))
    plt.plot(ySIS, yDeriv, '-', label=r"$\sigma$ = {}".format(p)) #label='di/dt~i'
    print("lamda={}\tmu={}\tsigma={}\t(1-1/sig)={}".format(lamda,lamda/p,p,(1-1/p)))

# 繪圖
plt.axhline(y=0,ls="--",c='c')  # 添加水平直線
plt.title("i(t)~di/dt in SIS model") # youcans-xupt
plt.legend(loc='best')
plt.show()

4. SIS 模型結果討論

SIS 模型代表：

若 $\sigma > 1$，則 $\lim\limits_{t \to \infty} i(t) = 1-1/\sigma$，代表患病者始終存在，成爲地方病。
若 $\sigma \leq 1$，則 $\lim\limits_{t \to \infty} i(t) = 0, (\sigma\leq 1)$ ，代表患病者人數不斷減小，最終能夠清零。
SIS 模型說明，對於傳染病，須要對患病者進行隔離以減小有效接觸，經過減小日接觸率 $\lambda$ 來減少接觸數 $\sigma$ ，打破傳播鏈，最終控制疫情。

須要指出的是，本文討論的 SIS模型是把考察地區視爲一個疫情均勻分佈的總體進行研究。實際上，在考察區域的疫情分佈必然是不均衡的，可能在局部區域發生疫情爆發致使該區域患病人數激增，是否會影響 SIS 模型的演化過程和穩定性呢？相關研究代表，擴散速度的不一樣可能致使種羣空間分佈的差別，在低風險區域將達到無病平衡點，在高風險區域仍將達到地方病平衡點。

【本節完】

歡迎關注『Python小白的數學建模課 @ Youcans』原創做品

原創做品，轉載必須標註原文連接：(https://www.cnblogs.com/youcans/p/14968504.html)。

Crated：2021-06-09

歡迎關注『Python小白的數學建模課 @ Youcans』，每週更新數模筆記
Python小白的數學建模課-01.新手必讀
 Python小白的數學建模課-02.數據導入
 Python小白的數學建模課-03.線性規劃
 Python小白的數學建模課-04.整數規劃
 Python小白的數學建模課-05.0-1規劃
 Python小白的數學建模課-06.固定費用問題
 Python小白的數學建模課-07.選址問題
 Python小白的數學建模課-09.微分方程模型
 Python小白的數學建模課-B2.新冠疫情 SI模型
 Python小白的數學建模課-B3.新冠疫情 SIS模型
 Python數模筆記-PuLP庫
 Python數模筆記-StatsModels統計迴歸
 Python數模筆記-Sklearn
Python數模筆記-NetworkX
Python數模筆記-模擬退火算法