大家還好嗎？

背景就不用多說了吧？本來我是初四上班的，現(xiàn)在延長(zhǎng)到2月10日了。這是我工作以來時(shí)間最長(zhǎng)的一個(gè)假期了?？上囊踩ゲ涣?。待在家里，沒啥事，就用python模擬預(yù)測(cè)一下新冠病毒肺炎的數(shù)據(jù)吧。要聲明的是本文純屬個(gè)人自娛自樂，不代表真實(shí)情況。

采用SIR模型，S代表易感者，I表示感染者，R表示恢復(fù)者。染病人群為傳染源，通過一定幾率把傳染病傳給易感人群，ta自己也有一定的幾率被治愈并免疫，或死亡。易感人群一旦感染即成為新的傳染源。

模型假設(shè):

①不考慮人口出生、死亡、流動(dòng)等情況，即人口數(shù)量保持常數(shù)。

②一個(gè)病人一旦與易感者接觸就必然具有一定的傳染力。假設(shè) t 時(shí)刻單位時(shí)間內(nèi)，一個(gè)病人能傳染的易感者數(shù)目與此環(huán)境內(nèi)易感者總數(shù)s(t)成正比，比例系數(shù)為β，從而在t時(shí)刻單位時(shí)間內(nèi)被所有病人傳染的人數(shù)為βs(t)i(t)。

③ t 時(shí)刻，單位時(shí)間內(nèi)從染病者中移出的人數(shù)與病人數(shù)量成正比，比例系數(shù)為γ，單位時(shí)間內(nèi)移出者的數(shù)量為γi(t)。
模型為

其中，β為感染系數(shù)，代表易感人群與傳染源接觸被感染的概率。γ為隔離(恢復(fù))系數(shù)，我們對(duì)其倒數(shù)1/γ更感興趣，代表了平均感染時(shí)間(average infectious period)。S(0)為初始易感人數(shù)，I(0)為初始感染人數(shù)。

按照[1]里面的代碼模型的感染人數(shù)是這樣的

現(xiàn)在的問題就是利用現(xiàn)有的數(shù)據(jù)找到新冠肺炎的β值，γ值等數(shù)據(jù)了。先把數(shù)據(jù)拔下來吧。從[3]上扒數(shù)據(jù)，由于數(shù)據(jù)不多，就手工完成吧。保存到csv文件里。

然后把數(shù)據(jù)作圖

還有一個(gè)指標(biāo)是再生數(shù)R0=β/γ，大于1時(shí)人群中大部分才被感染[4]。世衛(wèi)組織1月23日的估計(jì)是R0在1.4到2.5之間[5]，最新的根據(jù)前425例發(fā)病數(shù)據(jù)的估計(jì)值為2.2[6]。

文章[7]中的按一般病毒性肺炎恢復(fù)期25天計(jì)算得到的γ值為0.04。

關(guān)于β值和初始易感人群，[7]的作者采用的方法是先估計(jì)一個(gè)區(qū)間，然后用最小二乘法找到最佳參數(shù)，β≈3.57*10^-5。S[0]的范圍為5000-30000人。[7]文章里有matlab代碼，我用python改寫一下，由于對(duì)最小二乘法法的實(shí)現(xiàn)比較陌生，嘗試了半天，最后我決定用最笨的辦法——窮舉法。就是用兩個(gè)嵌套循環(huán)將范圍內(nèi)所有β值和S0值都試一遍，計(jì)算每次嘗試結(jié)果與實(shí)際數(shù)據(jù)之間差值的平方和，平方和最小的一組β值和S0值用來做預(yù)測(cè)。代碼如下:

γ值設(shè)定為0.04，即一般病程25天

用最小二乘法估計(jì)β值和初始易感人數(shù)

gamma = 0.04
S0 = [i for i in range(20000, 40000, 1000)]
beta = [f for f in np.arange(1e-7, 1e-4, 1e-7)]
# 定義偏差函數(shù)
def error(res):
 err = (data["感染者"] - res)**2
 errsum = sum(err)
 return errsum

# 窮舉法，找出與實(shí)際數(shù)據(jù)差的平方和最小的S0和beta值
minSum = 1e10
minS0 = 0.0
minBeta = 0.0
bestRes = None

for S in S0:
 for b in beta:
  # 模型的差分方程
  def diff_eqs_2(INP, t):
   Y = np.zeros((3))
   V = INP
   Y[0] = -b * V[0] * V[1]
   Y[1] = b * V[0] * V[1] - gamma * V[1]
   Y[2] = gamma * V[1]
   return Y

  # 數(shù)值解模型方程
  INPUT = [S, I0, 0.0]
  RES = spi.odeint(diff_eqs_2, INPUT, t_range)
  errsum = error(RES[:21, 1])
  if errsum < minSum:
   minSum = errsum
   minS0 = S
   minBeta = b
   bestRes = RES
   print("S0=%d beta=%f minErr=%f" % (S, b, errsum))
print("S0 = %d β = %f" % (minS0, minBeta))

結(jié)果 S0 = 39000, β = 8e-6

上述程序耗時(shí)較長(zhǎng)，只在探索時(shí)執(zhí)行，完了就注釋掉，用最優(yōu)參數(shù)進(jìn)行預(yù)測(cè)。

預(yù)測(cè)最大感染人數(shù):23769 時(shí)間是在1月10日的33天后，也就是2月12日。

本文代碼:https://github.com/zwdnet/2019-nCov-SIRmodel

與[7]作者討論，我的算法是將S0與β作為獨(dú)立的兩個(gè)變量用兩個(gè)循環(huán)嵌套分別遍歷，他的做法是用每個(gè)S0的值代入微分方程算出相應(yīng)的β值。他的算法應(yīng)該更好一些，我正在嘗試。另外在微信公眾號(hào)上看到一篇更系統(tǒng)的關(guān)于此次疫情的數(shù)學(xué)模型的文章:https://mp.weixin.qq.com/s/rgaJtA4jioLOCHs_oCauDg

再次聲明:本文只是我個(gè)人在家無聊的游戲作品，不是正兒八經(jīng)的預(yù)測(cè)。我也不是流行病學(xué)專業(yè)人士。祝疫情早日結(jié)束！武漢加油！中國(guó)加油！

總結(jié)

以上所述是小編給大家介紹的python模擬預(yù)測(cè)一下新型冠狀病毒肺炎的數(shù)據(jù)，希望對(duì)大家有所幫助！