Python小白的数学建模课-B3. 新冠疫情 SIS模型
傳染病的數學模型是數學建模中的典型問題,常見的傳染病模型有 SI、SIR、SIRS、SEIR 模型。
SIS 模型型將人群分為 S 類和 I 類,考慮患病者可以治愈而變成易感者,但不考慮免疫期。
本文詳細給出了 SIS 模型的建模、例程、運行結果和模型分析,讓小白都能懂。
『Python小白的數學建模課 @ Youcans』帶你從數模小白成為國賽達人。
Python小白的數學建模課-A3.12個新冠疫情數模競賽賽題及短評
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 模型適用于只有易感者和患病者兩類人群,且無法治愈的疾病。
按照 SI 模型,最終所有人都會被傳染而變成病人,這是因為模型中沒有考慮病人可以治愈。因此只能是健康人患病,而患病者不能恢復健康(甚至也不會死亡,而是不斷傳播疫情),所以終將全部被傳染。
1.2 SIS 模型
SIS 模型將人群分為 S 類和 I 類,考慮患病者(I 類)可以治愈而變成易感者(S 類),但不考慮免疫期,因此患病者(I 類)治愈變成易感者以后還可以被感染而變成患病者。
SIS 模型適用于只有易感者和患病者兩類人群,可以治愈,但會反復發作的疾病,例如腦炎、細菌性痢疾等治愈后也不具有免疫力的傳染病。
SIS 模型假設:
需要說明的是,不考慮生死或人口流動,通常是由于考慮一個封閉環境而且假定疫情隨時間的變化比生死、遷移隨時間的變化顯著得多, 因此后者可以忽略不計。
SIS 模型的微分方程:
由
Ndidt=Nλsi?NμiN\frac{di}{dt} = N \lambda s i - N \mu i Ndtdi?=Nλsi?Nμi
得:
didt=λi(1?i)?μi,i(0)=i0\frac{di}{dt} = \lambda i (1-i) - \mu i,\ i(0) = i_0 dtdi?=λi(1?i)?μi,?i(0)=i0?
由日治愈率 μ\muμ 可知平均治愈天數為 1/μ1/\mu1/μ,也稱平均傳染期。定義 σ=λ/μ\sigma = \lambda / \muσ=λ/μ,其含義是每個病人在傳染期內所傳染的平均人數,稱為傳染期接觸數。例如,平均傳染期 1/μ=51/\mu = 51/μ=5,日接觸率 λ=2\lambda = 2λ=2(每天傳染 2人),則傳染期接觸數 σ=10\sigma = 10σ=10。
SIS 模型的解析解為:
{i(t)=i01+λti0,λ=μi(t)=[λλ?μ+(1i0?λλ?μ)?e?(λ?μ)t]?1,λ≠μ\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}\\ ?????????i(t)=1+λti0?i0??i(t)=[λ?μλ?+(i0?1??λ?μλ?)?e?(λ?μ)t]?1?,λ=μ,λ?=μ??
注意:網上有些博文中解析解的公式誤寫成 exp((λ?μ)t)exp((\lambda-\mu)t)exp((λ?μ)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)f(y,t)f(y,t) ,即 y 在 t 處的導數,以函數的形式表示
- y0: array: 初始條件 y0y_0y0?,對于常微分方程組 y0y_0y0? 則為數組向量
- t: array: 求解函數值對應的時間點的序列。序列的第一個元素是與初始條件 y0y_0y0? 對應的初始時間 t0t_0t0?;時間序列必須是單調遞增或單調遞減的,允許重復值。
- args: 向導數函數 func 傳遞參數。當導數函數 f(y,t,p1,p2,..)f(y,t,p1,p2,..)f(y,t,p1,p2,..) 包括可變參數 p1,p2… 時,通過 args =(p1,p2,…) 可以將參數p1,p2… 傳遞給導數函數 func。
odeint() 的返回值:
- y: array 數組,形狀為 (len(t),len(y0),給出時間序列 t 中每個時刻的 y 值。
odeint() 的編程步驟:
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*ireturn 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 模型在 μ=0\mu = 0μ=0 時的特殊情況。
3. SIS 模型參數的影響
對于 SIS 模型,需要考慮日接觸率 λ\lambdaλ 與日治愈率 μ\muμ 的關系、患病者比例的初值 i0i_0i0? 的影響,總人數 N 沒有影響。
3.1 日接觸率 λ\lambdaλ 與日治愈率 μ\muμ 關系的影響
直觀地考慮,如果每天治愈的人數高于感染的人數,則疫情逐漸好轉,否則疫情逐漸嚴重。因此日接觸率 λ\lambdaλ 與日治愈率 μ\muμ 的關系非常關鍵,這就是傳染期接觸數 σ=λ/μ\sigma = \lambda / \muσ=λ/μ 的意義。
(1) σ≤1\sigma \leq 1σ≤1
當 σ<1\sigma<1σ<1 時,傳染期接觸數小于 1,日接觸率小于日治愈率,患病率單調下降,最終清零,與患病率初值無關。 σ\sigmaσ 越小,疫情清零速度越快; σ\sigmaσ 越接近于 1,疫情清零越慢,但最終仍將清零。
分析其實際意義,傳染期接觸數小于 1,表明在傳染期內經過接觸而使易感者變成患病者的數量,小于在傳染期內治愈的患病者的數量,因此患病者數量、比例都會逐漸降低,所以最終可以清零,稱為無病平衡點。
當 σ=1\sigma=1σ=1 時,不論患病率初值如何,患病率也是單調下降,最終趨近于 0。雖然在數學上患病率只能趨近于 0 而不等于 0,但考慮到總人數 N 是有限的,而患病者和易感者人數需要取整,因此 σ=1\sigma=1σ=1 時最終也會清零。
(2) σ>1\sigma > 1σ>1
當 σ>1\sigma>1σ>1 時,傳染期接觸數大于 1,日接觸率大于日治愈率,患病率的升降有兩種情況:
當患病率很低時,患病者人數少而易感者人數多,患病率上升;但隨著患病率增大,患病者越來越多而易感者越來越少,患病率雖然仍然上升但上升速度趨緩,最終趨于定值。
當患病率很高時,患病者人數多而易感者人數少,患病率下降;但隨著患病率減小,患病者越來越少而易感者越來越多,患病率雖然仍然下降但下降速度趨緩,最終也趨于相同的定值。
患病率最終都會收斂到穩態特征值 i∞=1?1/σi_\infty=1-1/\sigmai∞?=1?1/σ。當 i0>i∞i_0>i_\inftyi0?>i∞? 即患病率初值大于穩態特征值時,疫情曲線單調上升收斂;當 i0<i∞i_0<i_\inftyi0?<i∞? 即患病率初值小于穩態特征值時,疫情曲線單調下降收斂;當 i0=i∞i_0 = i_\inftyi0?=i∞? 時,患病率始終大于穩態特征值,疫情曲線為水平直線。
這表明,當 σ>1\sigma>1σ>1 時疫情終將穩定但不會清零,而是長期保持一定的患病率,稱為地方病平衡點。
當 σ=1\sigma=1σ=1 時,不論患病率初值如何,患病率都單調下降并最終趨于 0。
3.2 傳染期接觸數 σ\sigmaσ 與 di/dtdi/dtdi/dt 的關系
患病率的一階導數 di/dtdi/dtdi/dt 的變化曲線,表明不論傳染期接觸數和初值如何,患病率的變化率都將收斂到 0,因此疫情終將穩定。當 σ<1\sigma<1σ<1 時, di/dtdi/dtdi/dt 始終是負值,單調上升趨近于 0; 當 σ>1\sigma>1σ>1 時, di/dtdi/dtdi/dt 始終是正值,先上升達到峰值后再逐漸減小趨近于 0。
本圖為患病率 i(t)i(t)i(t) 與一階導數 di/dtdi/dtdi/dt 在不同傳染期接觸數下的關系曲線(相空間圖)。當 σ≤1\sigma\leq 1σ≤1 時,曲線收斂到原點 (0,0)(0,0)(0,0),即存在無病平衡點; 當 σ>1\sigma>1σ>1 時,曲線收斂到 (1?1/σ,0)(1-1/\sigma,0)(1?1/σ,0),即存在地方病平衡點。
3.3 Python例程:傳染期接觸數 σ\sigmaσ 與 di/dtdi/dtdi/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*ireturn 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 模型表明:
需要指出的是,本文討論的 SIS模型是把考察地區視為一個疫情均勻分布的整體進行研究。實際上,在考察區域的疫情分布必然是不均衡的,可能在局部區域發生疫情爆發導致該區域患病人數激增,是否會影響 SIS 模型的演化過程和穩定性呢?相關研究表明,擴散速度的不同可能導致種群空間分布的差異,在低風險區域將達到無病平衡點,在高風險區域仍將達到地方病平衡點。
【本節完】
版權聲明:
歡迎關注『Python小白的數學建模課 @ Youcans』 原創作品
原創作品,轉載必須標注原文鏈接:https://blog.csdn.net/youcans/article/details/117786272
Copyright 2021 Youcans, XUPT
Crated:2021-06-10
歡迎關注 『Python小白的數學建模課 @ Youcans』 系列,持續更新
Python小白的數學建模課-01.新手必讀
Python小白的數學建模課-02.數據導入
Python小白的數學建模課-03.線性規劃
Python小白的數學建模課-04.整數規劃
Python小白的數學建模課-05.0-1規劃
Python小白的數學建模課-06.固定費用問題
Python小白的數學建模課-07.選址問題
Python小白的數學建模課-09.微分方程模型
Python小白的數學建模課-10.微分方程邊值問題
Python小白的數學建模課-A1.國賽賽題類型分析
Python小白的數學建模課-A2.2021年數維杯C題探討
Python小白的數學建模課-A3.12個新冠疫情數模競賽賽題及短評
Python小白的數學建模課-B2. 新冠疫情 SI模型
Python小白的數學建模課-B3. 新冠疫情 SIS模型
Python小白的數學建模課-B4. 新冠疫情 SIR模型
Python小白的數學建模課-B5. 新冠疫情 SEIR模型
Python小白的數學建模課-B6. 新冠疫情 SEIR改進模型
Python數模筆記-PuLP庫
Python數模筆記-StatsModels統計回歸
Python數模筆記-Sklearn
Python數模筆記-NetworkX
Python數模筆記-模擬退火算法
總結
以上是生活随笔為你收集整理的Python小白的数学建模课-B3. 新冠疫情 SIS模型的全部內容,希望文章能夠幫你解決所遇到的問題。
- 上一篇: linux系统安装serv u,建立第一
- 下一篇: 【OpenCV 例程200篇】29. 图