Python實現新型冠狀病毒傳播模型及預測程式碼例項
1.傳染及發病過程
一個健康人感染病毒後進入潛伏期(時間長度為Q天),潛伏期之後進入發病期(時間長度為D天),發病期之後該患者有三個可能去向,分別是自愈、接收隔離、死亡。
2.模型假設
潛伏期Q=7天,根據報道潛伏期為2~14天,取中間值;發病期D=10天,根據文獻報告,WHO認定SARS發病期為10天,假設武漢肺炎與此相同;潛伏期的患者不具有將病毒傳染給他人的能力;發病期的患者具有將病毒傳染給他人的能力;患者在發病期之後不再具有將病毒傳染他人的能力;假設處於發病期的患者平均每天密切接觸1人,致使該人患病的概率為γ最初只有一個人類感染者;病情自然發展,沒有外部干擾。
3.模型公式
:人類感染該病毒的天數
: 第N天感染該病毒並且處於發病期的患者數量
: Q天前新被感染患者,當日進入發病期的數量
: 當日發病期滿,不再具有傳染能力的患者數量
4.模型初始值
根據假設,最初只有一個人類感染者,所以:
5.實際疫情資料
人類感染病毒且發病的初始日期:根據財新網的報道,官方通報首例不明原因肺炎是在12月8日,考慮到確診之前肯定已經嘗試過各種治療方案無效後認定為不明原因肺炎,所以有理由認為該名患者在12月8日已經處於發病期末端,根據假設發病期為10天,所以可以假設該名患者在11月29號發病,即N=1對應11月29日。
日新增發病數
近期疫情防控大事記:
從以上資訊可以判斷核酸檢測試劑是在1月16日、17日大幅使用的,18日、19日確診大量病例,因此1月20日之前的確診病例數對模型參考意義不大。1月20日之後,可以認為新發病例,發病即檢測。
截至1月20日24時,國家衛健委公告累計確診病例291
截至1月21日24時,國家衛健委公告累計確診病例440
截至1月22日24時,國家衛健委公告累計確診病例571
截至1月23日24時,國家衛健委公告累計確診病例830
截至1月24日24時,國家衛健委公告累計確診病例1287
所以:
1月21日新增確診病例:440-291=1491月22日新增確診病例:571-440=1311月23日新增確診病例:830-571=2591月24日新增確診病例:1287-830=457
考慮到:
1月20日之後發病即檢測確診檢測用時2日國外新增病例在個位數,且不能保證發病即檢測等因素,暫不考慮新增病例應該遞增,所以1月22日新增資料異常,捨去
使用上溯每日新增資料,同時考慮到確診需要2天, 可以得到:
,,
( 對應11月29日 )
6.擬合確定
根據近期每日新增資料、模型初始值及模型公式,用最小二乘擬合得到
7.預測
患者數量:根據上文確定的模型及引數,從11月29日(N=1)至1月27日(N=60)人群中累計處於發病期的人數如下圖所示:
根據模型,近期人群中患者數量計算如下:
人群中感染了病毒並處於發病期的患者數量,注意Pn一般需要延期2日才能確診
每日新增患者數量:
根據模型,近期人群中每日新增患者數量計算如下:
注意:圖中是人群中新增發病患者數量,可與晚2日的政府釋出新增數量進行對比。即22日新增患者數量可與24日政府釋出的新增病例進行對比。截至目前模型計算22日新增為369人,政府公佈的24日新增病例457人
根據以上模型預計未來幾日的情況如下:
1月25日將新確診432例,人群中發病患者為4068人;
1月26日將新確診505例,人群中發病患者為4759人;
1月27日將新確診590例,人群中發病患者為5568人;
1月28日將新確診691例,人群中發病患者為6514人;
1月29日將新確診809例,人群中發病患者為7621人。
由於1月20日之後採取的各種措施,將導致發病期D下降,感染概率 下降,1月29日之後日均增長勢頭會減弱。
8.模型程式碼
import numpy as np import matplotlib.pyplot as plt gamma = 0.55 Q = 7 D = 10 P = np.zeros(300,dtype=np.float) Psum = np.zeros(300,dtype=np.float) for i in range(Q): P[i] = 1 for j in range(300-Q): P[j+Q] = P[j+Q-1]+P[j]*gamma if j+Q-D-Q >= 0: P[j+Q] -= P[j+Q-D-Q]*gamma if j+Q == D: P[j+Q] -= 1 plt.xlabel("N") plt.ylabel("PN") plt.rcParams['font.sans-serif'] = ['SimHei'] plt.rcParams['axes.unicode_minus'] = False plt.plot(range(1,61),P[0:60]) plt.grid() plt.show()
以上就是我們小編給大家整理的全部內容,如果大家有任何補充可以聯絡小編。