跳轉到主要內容
工程的博客

貝葉斯建模使用PyMC3 COVID-19的時序動態

Srijith Rajamohan博士。
通過Srijith Rajamohan博士。

2021年1月6日 在工程的博客

分享這篇文章

在這篇文章中,我們看一下如何使用PyMC3為COVID-19推斷疾病的參數。PyMC3是一個流行的概率編程框架,用於貝葉斯建模。兩個受歡迎的方法來完成這是馬爾可夫鏈蒙特卡羅(密度),變分推理方法。工作在這裏看著感染病例使用目前可用的數據在美國作為一個時間序列並試圖模型使用一個區劃的概率模型。我們想嚐試推斷疾病最終估計參數R0使用獲得。

這裏介紹的工作是隻用於說明和現實生活中的貝葉斯建模需要更複雜的工具是什麼。關於種群動態是由各種假設,可能不是有效的大型非齊次的人群。同時,幹預措施,諸如社會距離,這裏不考慮接種疫苗。

這篇文章將介紹如下:

  1. 流行區劃的模型
  2. 數據從哪裏來,以及它是如何攝取嗎
  3. 先生,眾位模型的概述
  4. 貝葉斯推理與PyMC3常微分方程
  5. 推理工作流在磚

流行區劃的模型

區劃的模型和他們的行為的概述,請參閱這個筆記本在茱莉亞。

區劃的模型是一組常微分方程(ode)為封閉的數量,這意味著有一個運動的人口或隔間。這些目標人群的疾病傳播模型隔間均勻。你可以想象,這些假設可能不是有效的在龐大的人口。同樣重要的是要指出這至關重要的統計數據,如人口的出生和死亡人數可能不包含在這個模型。下麵提到的一些區劃的模型以及各種疾病傳播的隔間,然而,這不是一個詳盡的清單。

  • 容易感染(SI)
  • 容易感染恢複(先生)
  • 容易感染易感(SIS)
  • 容易感染恢複敏感(SIRS)
  • 容易感染恢複死亡(SIRD)
  • 容易受感染的恢複(西)
  • 容易受感染的易感恢複(西)
  • 容易暴露感染恢複死亡(SEIRD)
  • 母體免疫易感傳染性恢複(MSIR)
  • SIDARTHE

上麵列出的最後一個是最近的,專門針對COVID-19也許值得一讀的興趣。真實的疾病建模通常需要的不僅僅是疾病的演化階段由於許多與隔間違反相關的假設。了解疾病的傳播,我們想看的空間離散化和演化通過人口疾病的進展。框架的一個例子,這個時空演化模型是線(圖1)。

現實世界流行病建模(時空動力學)。

圖1

工具,如線使用人口普查數據和移動模式理解人們如何移動地理位置。線將全球劃分為空間網格大約25公裏的x 25公裏。大致有兩種類型的流動性:全球或遠程遷移和本地或短程遷移。長期流動主要包括航空旅行,因此機場被認為是疾病傳播的中心。海上旅行也是另一個重要因素,因此海軍港口是另一種類型的訪問點。隨著上述數學模型,這提供了一個隨機的框架,可以用來使數以百萬計的模擬得出推斷參數並作出預測。

數據來自約翰霍普金斯CSSE Github定期更新頁麵情況統計:

CSSE GitHub

例確診病例

的死亡人數

的數據可以從CSV文件可以通過Python讀入熊貓。

先生和先生們模型

爵士模型

爵士模型由一組三個常微分方程(ode)所示。在這個模型中有三個小櫃。

這裏‘S’,‘我’和‘R’指的是易感,感染和恢複部分的人口規模的“N”等

S +我+ R = N

這裏的假設是,一旦你已經從疾病中恢複過來,終生免疫是授予一個獨立的個體。這不是疾病的理由很多,因此可能不是一個有效的模型。

λ的感染率,μ的速度從疾病中康複。分數的人從感染中恢複是由“f”但是對於這項工作的目的,“f”設置為1。最後一個初值問題(IVP)為我們設置的常微分方程,我(0)被認為是已知的案例數量開始大流行和S(0)可以被估計為N -我(0)。這裏我們假設整個人口是敏感。我們的目標是做到以下幾點:

  • 利用貝葉斯推理估計λ和μ
  • 使用上述參數估計我任何時間(t) ' t '
  • 計算R0

正如已經指出的那樣,λ是透射係數的疾病。這取決於交互的數量,在單位時間內,感染人。這反過來又取決於感染的人口的數量。

λ=接觸率x傳輸概率

的力量感染或風險在任何時候“t”被定義為λΙt/Ν。另外,是經濟複蘇的一部分,在單位時間內發生的。μ1因此,平均恢複時間。基本的繁殖數量的R0是次要的平均數量情況下所產生的一個基本情況(例子嗎https://misuse.ncbi.nlm.nih.gov/error/abuse.shtml)。R0也定義的λ和μ的比率

R0=λ/μ(假設0接近1)

當R0> 1,我們有一個擴散的疾病大流行。最近努力接種脆弱,這已經變得更加相關的了解。如果我們“p”一小部分的人口接種疫苗(1 - p) R0眾位模型
眾位模型,如下所示,沒有假設的終生免疫感染者一旦恢複。因此,一個從恢複室隔間。因此,這可能是一個更好的低保真COVID-19基線模型,建議後天免疫是短期的。這裏唯一的附加參數γ是指免疫丟失的速度和受感染的個人從恢複池易感池。

對於這個工作,隻有先生模型實現,眾位模型和它的變體是留給未來的工作。

使用PyMC3推斷疾病的參數

先生我們可以離散化模型使用一階或二階時間分化方案,然後可以傳遞給PyMC3將3月時間提出的解決方案使用這些離散方程。參數λ和μ可以安裝使用蒙特卡羅抽樣程序。

一階方案

二級計劃

在PyMC3 DifferentialEquation方法

我們可以提供離散化時手動選擇一個高階離散化方案,這很快就會變得繁瑣且易出錯,更不用說計算效率低下。幸運的是,PyMC3有頌歌模塊幫助做的正是這一點。我們可以從ODE模塊使用DifferentialEquation方法將作為輸入一個函數返回的值的一組常微分方程是一個矢量,時間步驟所需的解決方案是,對應的狀態數的數量方程和變量的數量,我們想要解決。這種方法的缺點之一是,它往往是緩慢的。推薦的最佳實踐是使用“sunode”模塊在PyMC3(見下文)。例如,同樣的問題花費了5.4分鍾用DifferentialEquations與16 s sunode 100樣本,100年調樣品和20個時間點。

自我。sir_model_non_normalized = DifferentialEquation (func = self.SIR_non_normalized,* = self.time_range1:],n_states =2,n_theta =2,t0 =0)defSIR_non_normalized(自我,y, t, p):ds = - p [0]* y [0]* y [1)/ self.covid_data.N,di = p [0]* y [0]* y [1)/ self.covid_data。N - p (1]* y [1]返回(ds, di)

使用sunode模塊的語法如下所示。雖然有一些句法差異,DifferentialEquations的一般結構是一樣的。

進口sunode進口sunode.wrappers.as_theanodefSIR_sunode(t y p):返回{“年代”:- p。林* y。* y.I,“我”:p。林* y。* y。I - p。μ* y.I}
              ……
              sir_curves、_、問題解決者,_,_ = sunode.wrappers.as_theano.solve_ivp (y0 = {# ODE的初始條件“年代”:(S_init, ()),“我”:(I_init, ()),},params = {ODE的#參數,指定的形狀“林”(lam ()),“畝”(μ,()),“_dummy”:(np.array (1。),())#目前,sunode拋出一個錯誤},沒有這個號# RHS ODErhs = SIR_sunode,#時間點的解決方案tvals = *t0 = * [0),)

一個爵士模型的推理過程

為了對我們所尋求的參數進行推理,我們開始通過選擇合理的先驗參數的疾病。基於我們的理解這些參數的行為,對數正態分布是一個合理的之前。理想情況下,我們希望這個對數正態的平均參數的鄰居住我們預計所需的參數。良好的收斂性和解決方案,也是重要的數據可能是適當的(專業領域!)。通常選擇的可能性。

  • 正態分布
  • 對數正態分布
  • 學生的t分布

我們獲得敏感(S (t))和傳染性(我(t))從ODE求解器數字,然後樣本為λ和μ的值如下所示。

與pm.Model ()作為model4:σ= pm.HalfCauchy (“σ”,self.likelihood [“σ”),形狀=1)林= pm.Lognormal (“λ”,self.prior [“林”],self.prior [“lambda_std”])# 1.5,1.5μ= pm.Lognormal (“畝”,self.prior [“畝”],self.prior [“mu_std”])# 1.5,1.5res, _,問題解決者,_,_ = sunode.wrappers.as_theano.solve_ivp (y0 = {“年代”(自我。S_init ()),“我”(自我。I_init, ()),},params = {“林”(lam ()),“畝”(μ,()),“_dummy”:(np.array (1。),())},rhs = self.SIR_sunode,tvals = self.time_range,t0 = self.time_range [0])如果(可能性“分布”)= =對數正態的):我= pm.Lognormal (“我”μ= res (“我”σ=σ),觀察= self.cases_obs_scaled)elif(可能性“分布”)= =“正常”):我= pm.Normal (“我”μ= res (“我”σ=σ),觀察= self.cases_obs_scaled)elif(可能性“分布”)= =“students-t”):我= pm.StudentT (“我”ν=[可能性“怒”),#可能性分布的數據μ= res [“我”),#可能性分布的意思是,這些預測從先生σ=σ,觀察= self.cases_obs_scaled)R0 = pm.Deterministic (“R0”林/畝)跟蹤= pm.sample(自我。n_samples曲調=自我。n_tune鏈=4核=4)data = az.from_pymc3(跟蹤=跟蹤)

推理工作流與PyMC3磚

自開發一個這樣的模型,估算疾病參數使用貝葉斯推理,是一個迭代的過程我們希望盡可能自動化。這可能是一個好主意來實例化一個類模型對象的各種參數,自動運行。幸運的是,自動化執行使用磚筆記本很容易完成,每個單元格包含的組合所需的參數(見下文),一旦執行輸出塊無需用戶幹預。這也是一個好主意保存跟蹤信息,推理等指標R̂為每個運行以及其他元數據信息。文件格式如NetCDF可以用於這個雖然可以簡單使用Python內置數據庫模塊“擱置”。

covid_obj = COVID_data (“我們”、人口=328.2 e6)covid_obj.get_dates (data_begin =“10/1/20”data_end =“10/28/20”)sir_model = SIR_model_sunode (covid_obj)={可能性“分布”:“正常”,“σ”:2}之前= {“林”:1.5,“畝”:1.5,“lambda_std”:1.5,“mu_std”:1.5}sir_model.run_SIR_model (n_samples =500年n_tune =500年之前,可能性=可能性,=之前)

樣品結果

這些結果是純粹出於演示目的和廣泛的實驗之前,需要有意義的結果預計會從這個模擬。美國的病例數從1月到10月如下所示(圖2)。

例子COVID-19案例計算可視化生成PyMC3磚。

圖2

圖3顯示了一個推理的結果運行λ的後驗分布,μ和R0會顯示出來。執行貝葉斯推理的優點之一是,分布顯示中值估計隨著密度最高間隔(HDI)量化的不確定性。這是一個好主意來檢查跟蹤(至少),以確保抽樣是正確完成的。

例子一個推理運行結果顯示密度最高間隔使用PyMC3磚。

筆記和指導方針

一些建模和推理的一般規則:

  • 至少5000個樣本和1000個樣本用於調優
  • 上麵所示的結果,我已經使用:
    • 意思是:λ= 1.5μ= 1.5
    • 參數的標準偏差:2.0
  • 樣本至少3鏈
  • 集target_accept> 0.85
  • 如果可能的話,與示例核= n,“n”是可用的處理器數量
  • 檢查跟蹤收斂
  • 有限時域取樣影響推理的準確性,最好總是有更多的高質量的數據
  • 標準化數據,大值一般不收斂

調試你的模型

  • 自端PyMC3 theano, Python print語句不能用於檢查一個變量的值。使用theano.printing.Print (DESCRIPTIVE_STRING) (VAR)要做到這一點
  • 初始化隨機變量通過“testval”。這是非常有用的檢查那些討厭的壞能源的錯誤,通常是由於貧困的選擇可能或先知先覺。使用Model.check_test_point ()為了驗證這一點。
  • 使用一步= pm.Metropolis ()快速調試,運行更快但粗糙後的結果
  • 如果抽樣是緩慢的,之前檢查你的可能性分布

結論

這篇文章介紹了使用PyMC3獲得疾病的基本參數。在後續文章中,我們將看看如何使用磚環境和集成的工作流工具如MLflow實驗跟蹤和HyperOpt hyperparameter優化。

試著筆記本

引用

  • Priesemann集團工作
    • https://github.com/Priesemann-Group/covid_bayesian_mcmc
  • 德美特裏Pananos PyMC3頁麵
    • https://docs.pymc.io/notebooks/ODE_API_introduction.html
免費試著磚
開始
看到所有工程的博客的帖子
Baidu
map