這篇文章主要介紹python如何模擬預(yù)測一下新型冠狀病毒肺炎的數(shù)據(jù),文中介紹的非常詳細(xì),具有一定的參考價值,感興趣的小伙伴們一定要看完!
背景就不用多說了吧?本來我是初四上班的,現(xiàn)在延長到2月10日了。這是我工作以來時間最長的一個假期了??上囊踩ゲ涣?。待在家里,沒啥事,就用python模擬預(yù)測一下新冠病毒肺炎的數(shù)據(jù)吧。要聲明的是本文純屬個人自娛自樂,不代表真實情況。
采用SIR模型,S代表易感者,I表示感染者,R表示恢復(fù)者。染病人群為傳染源,通過一定幾率把傳染病傳給易感人群,ta自己也有一定的幾率被治愈并免疫,或死亡。易感人群一旦感染即成為新的傳染源。
模型假設(shè):
①不考慮人口出生、死亡、流動等情況,即人口數(shù)量保持常數(shù)。
②一個病人一旦與易感者接觸就必然具有一定的傳染力。假設(shè) t 時刻單位時間內(nèi),一個病人能傳染的易感者數(shù)目與此環(huán)境內(nèi)易感者總數(shù)s(t)成正比,比例系數(shù)為β,從而在t時刻單位時間內(nèi)被所有病人傳染的人數(shù)為βs(t)i(t)。
③ t 時刻,單位時間內(nèi)從染病者中移出的人數(shù)與病人數(shù)量成正比,比例系數(shù)為γ,單位時間內(nèi)移出者的數(shù)量為γi(t)。
模型為
其中,β為感染系數(shù),代表易感人群與傳染源接觸被感染的概率。γ為隔離(恢復(fù))系數(shù),我們對其倒數(shù)1/γ更感興趣,代表了平均感染時間(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ù)作圖
還有一個指標(biāo)是再生數(shù)R0=β/γ,大于1時人群中大部分才被感染[4]。世衛(wèi)組織1月23日的估計是R0在1.4到2.5之間[5],最新的根據(jù)前425例發(fā)病數(shù)據(jù)的估計值為2.2[6]。
文章[7]中的按一般病毒性肺炎恢復(fù)期25天計算得到的γ值為0.04。
關(guān)于β值和初始易感人群,[7]的作者采用的方法是先估計一個區(qū)間,然后用最小二乘法找到最佳參數(shù),β≈3.57*10^-5。S[0]的范圍為5000-30000人。[7]文章里有matlab代碼,我用python改寫一下,由于對最小二乘法法的實現(xiàn)比較陌生,嘗試了半天,最后我決定用最笨的辦法——窮舉法。就是用兩個嵌套循環(huán)將范圍內(nèi)所有β值和S0值都試一遍,計算每次嘗試結(jié)果與實際數(shù)據(jù)之間差值的平方和,平方和最小的一組β值和S0值用來做預(yù)測。代碼如下:
γ值設(shè)定為0.04,即一般病程25天
用最小二乘法估計β值和初始易感人數(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ù)據(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
上述程序耗時較長,只在探索時執(zhí)行,完了就注釋掉,用最優(yōu)參數(shù)進(jìn)行預(yù)測。
以上是“python如何模擬預(yù)測一下新型冠狀病毒肺炎的數(shù)據(jù)”這篇文章的所有內(nèi)容,感謝各位的閱讀!希望分享的內(nèi)容對大家有幫助,更多相關(guān)知識,歡迎關(guān)注創(chuàng)新互聯(lián)行業(yè)資訊頻道!
網(wǎng)站標(biāo)題:python如何模擬預(yù)測一下新型冠狀病毒肺炎的數(shù)據(jù)-創(chuàng)新互聯(lián)
鏈接URL:http://aaarwkj.com/article4/jddie.html
成都網(wǎng)站建設(shè)公司_創(chuàng)新互聯(lián),為您提供虛擬主機(jī)、微信小程序、品牌網(wǎng)站制作、網(wǎng)站營銷、靜態(tài)網(wǎng)站、網(wǎng)站排名
聲明:本網(wǎng)站發(fā)布的內(nèi)容(圖片、視頻和文字)以用戶投稿、用戶轉(zhuǎn)載內(nèi)容為主,如果涉及侵權(quán)請盡快告知,我們將會在第一時間刪除。文章觀點不代表本網(wǎng)站立場,如需處理請聯(lián)系客服。電話:028-86922220;郵箱:631063699@qq.com。內(nèi)容未經(jīng)允許不得轉(zhuǎn)載,或轉(zhuǎn)載時需注明來源: 創(chuàng)新互聯(lián)
猜你還喜歡下面的內(nèi)容