Past, Present, and Future of Software forBayesian Inference

貝葉斯推理軟件的過(guò)去、現(xiàn)在與未來(lái)

https://sites.stat.columbia.edu/gelman/research/published/Bayesian_software_review-8.pdf

摘要
在過(guò)去三十年中，隨著第一代 MCMC（馬爾可夫鏈蒙特卡洛）采樣器的普及，貝葉斯推斷的軟件工具經(jīng)歷了快速的發(fā)展。近年來(lái)，機(jī)器學(xué)習(xí)社區(qū)對(duì)新軟件包的積極開發(fā)以及某些特定應(yīng)用領(lǐng)域?qū)Ｓ密浖牧餍校M(jìn)一步刺激了用戶數(shù)量的指數(shù)增長(zhǎng)。本綜述旨在總結(jié)目前最流行的貝葉斯推斷軟件，并為讀者提供一份實(shí)用指南，幫助其在貝葉斯計(jì)算的世界中導(dǎo)航。我們預(yù)計(jì)，在諸如概率編程、無(wú)需似然的推斷以及貝葉斯神經(jīng)網(wǎng)絡(luò)等多個(gè)研究領(lǐng)域中，相關(guān)算法及其配套軟件將繼續(xù)快速發(fā)展，這將進(jìn)一步拓展貝葉斯范式在各類激動(dòng)人心的應(yīng)用中的可能性。

關(guān)鍵詞和短語(yǔ) ：統(tǒng)計(jì)學(xué)、數(shù)據(jù)分析、MCMC、計(jì)算、概率編程

1. 引言

在過(guò)去三十年中，貝葉斯推斷已經(jīng)確立了其作為經(jīng)典統(tǒng)計(jì)推斷方法的一種可行替代方案的地位，并如今已成為每一位統(tǒng)計(jì)學(xué)家工具箱中必備的工具。許多理論和方法上的發(fā)展都促成了貝葉斯統(tǒng)計(jì)的成功。然而，沒有任何一項(xiàng)進(jìn)展對(duì)于大眾化采用貝葉斯方法的重要性能夠超過(guò)那些易于使用且穩(wěn)健的軟件 的出現(xiàn)。

我們撰寫本文的目標(biāo)是向讀者介紹貝葉斯推斷軟件的歷史、現(xiàn)狀與未來(lái)發(fā)展趨勢(shì) 。我們希望為讀者提供一份全面的綜述，涵蓋流行的軟件工具、推動(dòng)這些軟件發(fā)展的統(tǒng)計(jì)學(xué)與計(jì)算技術(shù)的關(guān)鍵進(jìn)展，以及當(dāng)前所面臨的局限性和挑戰(zhàn)。本文既面向從事貝葉斯統(tǒng)計(jì)實(shí)踐的專業(yè)人士，也面向?qū)@一領(lǐng)域不太熟悉、但希望了解貝葉斯推斷任務(wù)及其所用工具的讀者。

在正式展開之前，我們先簡(jiǎn)要回顧一下背景知識(shí)，并介紹一些在全文中將頻繁使用的術(shù)語(yǔ)。

1.1 貝葉斯推斷

貝葉斯推斷的核心思想是：將選定的似然函數(shù) p(y∣θ)和參數(shù) θ的先驗(yàn)分布 p(θ)（或聯(lián)合分布 p(θ,y)）與觀測(cè)數(shù)據(jù) y結(jié)合起來(lái)，從而計(jì)算出參數(shù)的后驗(yàn)分布 p(θ∣y)。我們通過(guò)貝葉斯定理 來(lái)完成這一過(guò)程：

在貝葉斯推斷中，最常見的關(guān)注量是參數(shù)或其函數(shù)的 后驗(yàn)性質(zhì) ，這些性質(zhì)可以通過(guò)對(duì)后驗(yàn)分布 p(θ∣y) 求期望的方式來(lái)表達(dá)：

先驗(yàn)預(yù)測(cè)和后驗(yàn)預(yù)測(cè)、模型選擇以及其他關(guān)注量遵循類似的模式。因此，貝葉斯推斷中的主要計(jì)算問(wèn)題就是計(jì)算積分 。

我們所選擇的似然函數(shù)和先驗(yàn)分布很少能為 p(θ∣y)提供閉合形式的解，在大多數(shù)情況下，我們只能將它計(jì)算到一個(gè)乘法常數(shù)的精度；而對(duì)于積分本身，就更少有閉合解了。因此，計(jì)算這些關(guān)注量本質(zhì)上是一個(gè)數(shù)值問(wèn)題，本身就具有挑戰(zhàn)性。

1.2 貝葉斯推斷軟件

在討論貝葉斯推斷軟件時(shí)，我們將軟件組件分為三類：建模語(yǔ)言 （modeling language）、計(jì)算方法 （computation methods） 和 實(shí)用工具 （utilities）。

1.2.1 建模語(yǔ)言

我們以最廣義的方式來(lái)使用“建模語(yǔ)言”這一術(shù)語(yǔ)，指代一種允許用戶指定似然函數(shù)、先驗(yàn)分布和數(shù)據(jù)的組件（從現(xiàn)在起，我們用“模型”一詞來(lái)統(tǒng)稱所有這些內(nèi)容）。另一種方式是通過(guò)指定聯(lián)合分布 p(θ,y)的生成模型來(lái)進(jìn)行貝葉斯推斷（見第4.3節(jié)），有些語(yǔ)言支持同時(shí)指定這兩種方式。附錄中給出了不同建模語(yǔ)言的一個(gè)示例說(shuō)明。

每種建模語(yǔ)言都在通用性與易用性 之間做出某種權(quán)衡。在光譜的一端是一些表達(dá)能力很強(qiáng)的語(yǔ)言，例如通用的概率編程語(yǔ)言（PPLs）以及像 Python 這樣的通用編程語(yǔ)言；而在另一端，則是一些僅支持單一模型或有限選項(xiàng)的軟件；介于兩者之間的則是專門為貝葉斯推斷設(shè)計(jì)的聲明式語(yǔ)言（如 WinBUGS [74]）、命令式語(yǔ)言（如 Stan [21]）或基于公式的語(yǔ)言（如 R-INLA [71] 和 rstanarm [6]），后者使用的語(yǔ)法類似于 R 核心統(tǒng)計(jì)包 stats 中廣義線性模型（GLMs）所使用的公式對(duì)象語(yǔ)法 [121] 等。

與其他任何組件相比，建模語(yǔ)言的選擇 更能決定目標(biāo)用戶群體。或者反過(guò)來(lái)說(shuō)，當(dāng)軟件是為特定用戶群體設(shè)計(jì)時(shí)，沒有任何其他組件會(huì)像建模語(yǔ)言那樣受到目標(biāo)用戶需求的深刻影響。正如各種不同的建模語(yǔ)言所展示的那樣，貝葉斯推斷的用戶群體是異質(zhì)性的，并不存在一種適用于所有人的統(tǒng)一方案。

1.2.2 計(jì)算方法

一旦模型被指定，下一步就是執(zhí)行對(duì)后驗(yàn)分布以及其他關(guān)注量的計(jì)算。因此，一個(gè)完整的貝葉斯推斷軟件必須實(shí)現(xiàn)一種或多種貝葉斯計(jì)算方法。

沒有任何一種方法可以為所有模型都實(shí)現(xiàn)實(shí)際可行的貝葉斯計(jì)算。因此，人們開發(fā)了多種不同的計(jì)算方法，每種方法都在通用性與效率之間進(jìn)行權(quán)衡 。計(jì)算方法決定了哪些類別的模型是可以被計(jì)算的，并且通常比建模語(yǔ)言更限制軟件的能力。也就是說(shuō)，建模語(yǔ)言可能允許用戶定義某些模型，但這些模型在理論上甚至都無(wú)法被計(jì)算方法處理的情況并不少見。而且一般來(lái)說(shuō)，總存在一些模型，在實(shí)踐中某個(gè)計(jì)算方法無(wú)法處理，即使它在理論上是可行的。

在本文中，我們從貝葉斯軟件的角度 來(lái)討論貝葉斯計(jì)算：我們只限于討論那些對(duì)貝葉斯推斷軟件發(fā)展至關(guān)重要的方法，并列出各個(gè)軟件中所實(shí)現(xiàn)的方法。關(guān)于貝葉斯計(jì)算的歷史和最新進(jìn)展的詳細(xì)內(nèi)容，我們推薦閱讀文獻(xiàn) [81]。

1.2.3 實(shí)用工具（Utilities）

我們將“實(shí)用工具”定義為不屬于前兩類（建模語(yǔ)言和計(jì)算方法）的所有軟件組件。盡管如此，它們?cè)谪惾~斯軟件中非常常見，并且對(duì)于貝葉斯推斷的工作流程來(lái)說(shuō)，即便不是必不可少的，也是十分便利的（關(guān)于貝葉斯工作流程的詳細(xì)討論，可參閱 [41, 46]）：

• 貝葉斯計(jì)算診斷 ：貝葉斯計(jì)算方法可能會(huì)失敗，比如無(wú)法找到最優(yōu)解；或者在使用馬爾可夫鏈蒙特卡洛（MCMC）方法時(shí)，無(wú)法正確地探索后驗(yàn)分布。診斷工具對(duì)于在解釋結(jié)果之前識(shí)別潛在問(wèn)題至關(guān)重要。此外，大多數(shù)關(guān)鍵方法都是基于 MCMC 的，因此是基于抽樣且近似的。近似誤差也必須被量化并在結(jié)果解釋中加以考慮。常見的診斷方法包括追蹤圖（traceplots）、蒙特卡洛標(biāo)準(zhǔn)誤、有效樣本量（ESS）、R? 統(tǒng)計(jì)量以及基于仿真的校準(zhǔn)（simulation-based calibration）[115]。

• 模型驗(yàn)證與比較 ：包括先驗(yàn)與后驗(yàn)可視化、先驗(yàn)預(yù)測(cè)檢查與后驗(yàn)預(yù)測(cè)檢查、（近似）留一交叉驗(yàn)證（leave-one-out cross-validation）以及模型評(píng)估準(zhǔn)則如 WAIC [134]，還有貝葉斯因子（Bayes factors）的計(jì)算等。建模語(yǔ)言決定了這些任務(wù)的難易程度。例如，為了進(jìn)行先驗(yàn)和后驗(yàn)預(yù)測(cè)檢查，我們需要從先驗(yàn)分布 p(y)和后驗(yàn)預(yù)測(cè)分布中抽樣；為了計(jì)算貝葉斯因子，我們需要評(píng)估邊緣似然 p(y)；而為了進(jìn)行交叉驗(yàn)證，則需要評(píng)估后驗(yàn)預(yù)測(cè)分布 。
  

• 計(jì)算庫(kù) ：包括矩陣代數(shù)庫(kù)、概率分布和其他統(tǒng)計(jì)計(jì)算的支持、高性能計(jì)算支持，以及自動(dòng)微分（AD）庫(kù)。

• 接口 ：通常，貝葉斯軟件只向用戶提供底層的命令行界面，用戶需要將數(shù)據(jù)和模型以文件形式傳入。為了使用方便，人們開發(fā)了更高級(jí)的接口，使用戶能夠通過(guò)流行的高級(jí)語(yǔ)言（如 Python 和 R）來(lái)調(diào)用這些計(jì)算功能。

• 文檔 ：包括軟件文檔、語(yǔ)言定義、示例、案例研究以及其他有助于使用該軟件的資料。

1.3 范圍與結(jié)構(gòu)安排

本文部分地綜述了目前最流行、以及歷史上最重要的通用貝葉斯推斷軟件。我們也包括了一些具有更特定用途的流行軟件。例如，僅提供某種實(shí)用功能的貝葉斯計(jì)算工具，或?qū)Ｗ⒂谀骋活愑邢弈Ｐ偷能浖?/p>

對(duì)于那些使用較少且用途更具體的軟件，本文的介紹偏向于 Python 和 R 這兩種最流行的數(shù)據(jù)分析語(yǔ)言。請(qǐng)參見表1和表2，分別對(duì) Python 和 R 中貝葉斯軟件包的相對(duì)流行程度進(jìn)行了估計(jì)。

通用貝葉斯計(jì)算經(jīng)歷了兩個(gè)明顯不同的階段，每個(gè)階段都由某一種類型的貝葉斯計(jì)算方法及其配套軟件主導(dǎo)：

• 從1990年代初到2010年代

  ，主導(dǎo)方法是 Gibbs 采樣 ，其代表性軟件是 BUGS ；

• 從2010年代至今

  ，主導(dǎo)方法是 哈密頓蒙特卡洛（HMC） ，其代表性軟件是 Stan 。

本文余下部分的第一部分大致對(duì)應(yīng)這兩個(gè)時(shí)期：

• 第2節(jié) 描述了 Gibbs 采樣、基于 Gibbs 的軟件典型結(jié)構(gòu)，以及 BUGS 語(yǔ)言。我們還介紹了可能在之后開發(fā)但與 BUGS 相關(guān)、受其啟發(fā)或作為其延續(xù)的軟件。
• 第3節(jié) 則聚焦于 HMC 和 Stan。

我們將 第4節(jié) 專門用于那些無(wú)法明確歸入上述兩個(gè)時(shí)期的軟件。這部分涵蓋了以計(jì)算為核心、面向特定模型類別，以及貝葉斯軟件和通用概率編程語(yǔ)言（PPLs）最新發(fā)展的工具。

最后，在 第5節(jié) 中，我們討論了貝葉斯推斷軟件的未來(lái)發(fā)展方向。

2. 第一代 —— 基于 Gibbs 采樣的方法

在20世紀(jì)90年代初到21世紀(jì)初的這段時(shí)間里，最流行的通用貝葉斯推斷軟件是基于圖模型 和Gibbs 采樣 作為貝葉斯計(jì)算方法的。

這種方法的主要假設(shè)是：我們聯(lián)合分布 p(V)=p(θ,y)中變量之間的條件獨(dú)立性 可以用一個(gè)有向無(wú)環(huán)圖（DAG） 來(lái)表示。在這個(gè)圖中，每個(gè)變量由一個(gè)節(jié)點(diǎn)表示，并且在給定其馬爾可夫毯（Markov blanket）的條件下，該節(jié)點(diǎn)與其他所有節(jié)點(diǎn)條件獨(dú)立。

能夠用這種圖結(jié)構(gòu)表示的模型被稱為貝葉斯網(wǎng)絡(luò) （Bayesian network），它是概率圖模型（probabilistic graphical models）的一個(gè)類別（附錄中提供了一個(gè)示例）。此時(shí)，聯(lián)合分布可以分解為如下形式：

一個(gè)每次更新一個(gè)節(jié)點(diǎn) 、并使用其完整條件分布 （full conditional）的馬爾可夫鏈，在較弱的條件下會(huì)收斂到后驗(yàn)分布。從實(shí)際角度來(lái)看，這意味著我們只需要能夠迭代地從完整條件分布中抽樣 即可。

算法1是對(duì) Gibbs 采樣算法的一個(gè)總結(jié)。該算法的一大優(yōu)勢(shì)在于不需要調(diào)整任何算法參數(shù) ，這對(duì)于自動(dòng)化推斷來(lái)說(shuō)是一個(gè)非常有用的特性。

在從完整條件分布中進(jìn)行抽樣的過(guò)程中，通常會(huì)使用一個(gè)抽樣器的層次結(jié)構(gòu) （hierarchy of samplers）。由于模型是以符號(hào)形式表達(dá)的，因此很容易檢查完整條件分布的性質(zhì)。在大多數(shù)情況下，分布的形式越具體，我們可以使用的抽樣算法就越高效 。

2.1 BUGS

這一方法的代表性軟件是 BUGS （Bayesian inference Using Gibbs Sampling，即“使用 Gibbs 采樣的貝葉斯推斷”）[110]。BUGS 項(xiàng)目于 1989 年在劍橋大學(xué)醫(yī)學(xué)研究委員會(huì)生物統(tǒng)計(jì)學(xué)系啟動(dòng)。BUGS 軟件后來(lái)演變?yōu)?WinBUGS [74, 76, 112]，它更新了 BUGS 語(yǔ)言和采樣算法；以及 OpenBUGS ，它是 WinBUGS 的 GNU 通用公共許可證（GPL）版本，并可在 Linux 上運(yùn)行（有一定限制）[111]。

目前，BUGS、WinBUGS 和 OpenBUGS 已不再積極開發(fā)，但 BUGS 項(xiàng)目啟發(fā)了許多其他軟件，我們將在本節(jié)末尾討論這些后續(xù)發(fā)展。BUGS 的詳細(xì)歷史可參見 Lunn 等人 [75]。

公式 (1) 所示的聯(lián)合分布分解形式，既是 BUGS 語(yǔ)言的基礎(chǔ) ，也是基于 Gibbs 采樣的計(jì)算方法的基礎(chǔ)。BUGS 語(yǔ)言是一種聲明式語(yǔ)言 ，用戶在其中聲明模型中變量之間的所有父子關(guān)系 P(v)。附錄中給出了一個(gè)用 JAGS 編寫的貝葉斯網(wǎng)絡(luò)模型示例，JAGS 是一種與 WinBUGS 非常相似的語(yǔ)言。

為了從完整條件分布中進(jìn)行抽樣，WinBUGS 針對(duì)以下不同情況分別實(shí)現(xiàn)了不同的方法 [74, 第12.1章]：

• 離散分布 ：使用逆累積分布函數(shù)法（inverse CDF method）；
• 標(biāo)準(zhǔn)分布 ：使用針對(duì)該分布的標(biāo)準(zhǔn)算法；
• 對(duì)數(shù)凹分布（log-concave） ：使用無(wú)需導(dǎo)數(shù)的自適應(yīng)拒絕抽樣（adaptive rejection sampling）[47]。許多標(biāo)準(zhǔn)分布都是對(duì)數(shù)凹的，包括指數(shù)族分布。多個(gè)對(duì)數(shù)凹分布的乘積仍然是對(duì)數(shù)凹的，因此完整條件分布通常也具有對(duì)數(shù)凹性質(zhì)；
• 受限范圍分布 ：使用切片抽樣（slice sampling）[91]；
• 無(wú)約束范圍分布 ：使用當(dāng)前點(diǎn) Metropolis 算法。

OpenBUGS 還引入了塊抽樣方法 （block sampling），可以同時(shí)對(duì)根據(jù)模型結(jié)構(gòu)判斷可能相關(guān)的一組節(jié)點(diǎn)進(jìn)行聯(lián)合抽樣。塊更新解決了 Gibbs 抽樣的一個(gè)主要缺點(diǎn)：其性能強(qiáng)烈依賴于模型的參數(shù)化方式。如果兩個(gè)變量后驗(yàn)相關(guān)性很高，但在 Gibbs 抽樣中被獨(dú)立更新，那么馬爾可夫鏈對(duì)于這兩個(gè)變量都會(huì)表現(xiàn)出高度自相關(guān)。通過(guò)對(duì)相關(guān)節(jié)點(diǎn)進(jìn)行塊更新可以解決這個(gè)問(wèn)題，否則用戶只能通過(guò)重新參數(shù)化模型來(lái)應(yīng)對(duì)這一問(wèn)題。

BUGS 概率編程語(yǔ)言的一個(gè)強(qiáng)項(xiàng)在于：數(shù)據(jù)與參數(shù)的區(qū)別是在運(yùn)行時(shí)根據(jù)提供的觀測(cè)值自動(dòng)識(shí)別的 。向量也可以部分觀測(cè)，未觀測(cè)的元素只需設(shè)為未知（NA）。這簡(jiǎn)化了用于后驗(yàn)檢查的抽樣模擬過(guò)程。

雖然 WinBUGS 主要面向貝葉斯網(wǎng)絡(luò)，但它也對(duì)無(wú)向圖 （如因子模型）提供了一定支持，前提是將變量的整個(gè)子集表示為一個(gè)單獨(dú)的多元節(jié)點(diǎn)，從而對(duì)其值進(jìn)行聯(lián)合抽樣。WinBUGS 還通過(guò) DoodleBUGS 編輯器支持使用 plate notation（板式表示法）進(jìn)行圖形化模型定義。

MultiBUGS [54] 是 BUGS 項(xiàng)目的延續(xù)。它的主要貢獻(xiàn)在于提供了更高效的實(shí)現(xiàn)方式以及多種并行化技術(shù)。在多核環(huán)境中，MultiBUGS 的效率可以比 OpenBUGS 高出幾個(gè)數(shù)量級(jí)。

此外，還有針對(duì) WinBUGS、OpenBUGS 和 MultiBUGS 的 R 接口：

• R2WinBUGS [113]
• R2OpenBUGS [126]
• R2MultiBUGS

JAGS（Just Another Gibbs Sampler） [97] 在語(yǔ)言和計(jì)算方式上與 WinBUGS 類似（差異見 [74, 第12.6章]）。不同于用 Component Pascal 編寫的 WinBUGS 和 OpenBUGS，JAGS 是用 C++ 編寫 的，并具有良好的可移植性。這促進(jìn)了它的流行，并且它目前仍在積極開發(fā)中。

附錄中提供了一個(gè)用 JAGS 編寫的模型示例。

JAGS 是 BUGS 的一個(gè)“克隆版”，擁有完全獨(dú)立的代碼庫(kù)，但目標(biāo)是實(shí)現(xiàn)類似的功能（詳見 [74, 第12.6章] 的差異總結(jié)）。JAGS 采用 C++ 編寫，支持 Windows、MacOS 和 Linux 系統(tǒng)，發(fā)布于 GNU 通用公共許可證第2版之下。

JAGS 內(nèi)嵌了一份 R 數(shù)學(xué)庫(kù)的副本，提供了高質(zhì)量的隨機(jī)數(shù)生成算法以及與概率分布相關(guān)的各種計(jì)算功能。JAGS 的主要采樣方法是 切片抽樣（slice sampling） [91]，該方法適用于連續(xù)變量和離散變量節(jié)點(diǎn)。

JAGS 的 “glm” 模塊包含針對(duì)廣義線性混合模型 （GLMMs）的高效采樣器。這些采樣器基于數(shù)據(jù)增強(qiáng) （data augmentation）技術(shù)，這是一種常用的通過(guò)添加新節(jié)點(diǎn)來(lái)簡(jiǎn)化圖模型抽樣的技巧 [59]。在這種情況下，數(shù)據(jù)增強(qiáng)將具有二元結(jié)果的 GLMM [1, 57, 98] 或同時(shí)具有二元和泊松結(jié)果的 GLMM [39] 轉(zhuǎn)換為具有正態(tài)結(jié)果的線性模型。這種轉(zhuǎn)換使得可以對(duì)線性預(yù)測(cè)器中的所有參數(shù)進(jìn)行塊更新 （block updating），這比傳統(tǒng)的 Gibbs 抽樣要高效得多。

該線性模型的底層引擎使用了稀疏矩陣代數(shù) [29]，能夠同時(shí)處理固定效應(yīng)和隨機(jī)效應(yīng)。

2.3 Nimble

Nimble [30] 與 BUGS 類似，也專注于圖模型。它是 BUGS 語(yǔ)言的一個(gè)擴(kuò)展，同時(shí)還實(shí)現(xiàn)了一種嵌入在 R 中的建模語(yǔ)言 ，兩者都會(huì)被編譯為 C++ 代碼。

Nimble 實(shí)現(xiàn)了多種貝葉斯計(jì)算方法，包括：

• Metropolis-Hastings（MH）

• Gibbs 抽樣

• 序貫蒙特卡洛（Sequential Monte Carlo，SMC）

用戶具有靈活性，可以為不同的節(jié)點(diǎn)分配不同的采樣方法。最近，Nimble 還增加了對(duì)自動(dòng)微分（AD）和哈密頓蒙特卡洛（HMC）的支持。

附錄中提供了一個(gè)用 Nimble 的 R 嵌入式語(yǔ)言編寫的模型示例。

為了模擬哈密頓動(dòng)力學(xué)，我們需要以某個(gè)步長(zhǎng) ?對(duì)時(shí)間進(jìn)行離散化。最常用的方法是蛙跳辛積分器（leapfrog symplectic integrator） 。

哈密頓動(dòng)力學(xué)具有幾個(gè)對(duì)于 HMC 能夠正常運(yùn)行至關(guān)重要的性質(zhì)：

• 它保持哈密頓量 H 不變，

• 它是可逆的（reversible），

• 它是辛的（symplectic），因此也保持體積不變（volume preserving）。

在 HMC 中，通常選擇哈密頓量 H使其具有可分離性（separable）形式：

H(q,p)=U(q)+K(p)

其中 U(q)是系統(tǒng)的勢(shì)能，K(p)是動(dòng)能。

HMC 的核心思想是利用哈密頓量來(lái)定義位置和動(dòng)量的聯(lián)合密度：

聯(lián)合密度 p(q,p)可以看作是在位置向量 q的目標(biāo)密度基礎(chǔ)上，通過(guò)為動(dòng)量向量 p添加一個(gè)獨(dú)立的多元高斯分布（均值為 0，協(xié)方差矩陣為 M）而得到的增廣密度。

哈密頓動(dòng)力學(xué)保持哈密頓量不變，因此軌跡上的所有狀態(tài)都具有相同的密度 p(?,?)。這使得哈密頓動(dòng)力學(xué)非常適合用于在 MCMC 算法中提出下一個(gè)狀態(tài)，因?yàn)橐粭l軌跡可以提出一個(gè)在位置 q上距離當(dāng)前狀態(tài)很遠(yuǎn)的新狀態(tài)，但接受概率仍然為 1。

為了能夠訪問(wèn)到所有可能的狀態(tài)，我們必須每次采樣一個(gè)新的動(dòng)量。由于聯(lián)合密度中的動(dòng)能和勢(shì)能部分是相互獨(dú)立的，并且我們是從動(dòng)量 p的真實(shí)分布中進(jìn)行采樣的，這種采樣不會(huì)改變目標(biāo)分布。也就是說(shuō)，p(q,p)仍然是馬爾可夫鏈的平穩(wěn)分布。

然而在實(shí)際中，盡管蛙跳積分方法對(duì)哈密頓動(dòng)力學(xué)的模擬是穩(wěn)定的，但它并不能完全保持哈密頓量不變，會(huì)出現(xiàn)相對(duì)較小的波動(dòng)。這就是為什么我們?nèi)匀恍枰獞?yīng)用Metropolis 校正 的原因。

將上述步驟綜合起來(lái)，算法2總結(jié)了基本的 HMC 算法。

HMC 背后的核心思想早在其被廣泛應(yīng)用于貝葉斯軟件之前就已經(jīng)為人所知超過(guò) 20 年 [33]。使 HMC 更加自動(dòng)化使用的關(guān)鍵推動(dòng)力是自動(dòng)微分 （AD，automatic differentiation）的發(fā)展（參見 [5] 對(duì)自動(dòng)微分的介紹與綜述）。

模擬 HMC 的哈密頓動(dòng)力學(xué)需要目標(biāo)密度的梯度信息。為了使軟件具備通用性，我們必須能夠?yàn)橛脩艟帉懙乃谐绦蛴?jì)算出梯度。在四種常見的導(dǎo)數(shù)計(jì)算方法中，有三種并不可行：

• 手動(dòng)推導(dǎo)導(dǎo)數(shù)不具實(shí)用性；

• 通過(guò)有限差分進(jìn)行數(shù)值微分由于舍入誤差和截?cái)嗾`差而過(guò)于不穩(wěn)定，并且在高維情況下效率低下；

• 符號(hào)微分則會(huì)遇到表達(dá)式膨脹問(wèn)題（expression swell），導(dǎo)致生成的代碼效率低下。

而自動(dòng)微分（AD）利用了這樣一個(gè)事實(shí)：每個(gè)程序都是由基本運(yùn)算組成的；只要每一個(gè)基本運(yùn)算都實(shí)現(xiàn)了其導(dǎo)數(shù)，我們就可以應(yīng)用鏈?zhǔn)椒▌t來(lái)求得整個(gè)程序的梯度。這使得我們可以以機(jī)器精度 計(jì)算梯度。現(xiàn)代大多數(shù)推斷軟件都內(nèi)置或引入了一個(gè)自動(dòng)微分庫(kù)。

HMC 的一個(gè)重要局限性在于：它只能用于連續(xù)空間 （smooth spaces）。

要讓 HMC 在通用貝葉斯推斷中發(fā)揮作用，一個(gè)關(guān)鍵挑戰(zhàn)是自動(dòng)調(diào)參 ，包括質(zhì)量矩陣（mass matrix）、步長(zhǎng)（step size）以及步數(shù)（number of steps）。基于 HMC 的軟件通常會(huì)實(shí)現(xiàn)一個(gè)或多個(gè)預(yù)熱階段（warmup phase）來(lái)進(jìn)行參數(shù)調(diào)優(yōu)。之后進(jìn)入采樣階段，預(yù)熱期間的樣本會(huì)被丟棄。

其中一項(xiàng)關(guān)鍵發(fā)展是 NUTS（無(wú) U 型轉(zhuǎn)彎采樣器） [61]，它經(jīng)過(guò)一些改進(jìn)后仍然是 Stan 中的核心貝葉斯計(jì)算方法。NUTS 的基本思想是動(dòng)態(tài)調(diào)整步數(shù)：持續(xù)模擬哈密頓軌跡，直到檢測(cè)到軌跡開始“掉頭”返回初始狀態(tài)（或達(dá)到最大步數(shù)限制）為止。

盡管存在一些有前景的替代方案可用于調(diào)節(jié)步數(shù)（包括更適合 GPU 計(jì)算的變體版本）[60, 62]，但 NUTS 仍是目前最常用的 HMC 實(shí)現(xiàn)方式，并且也被集成到了大多數(shù)現(xiàn)代通用貝葉斯推斷軟件中。

HMC / NUTS 支持幾種特定的 MCMC 診斷工具 [10]：

• 當(dāng)步長(zhǎng)過(guò)大，無(wú)法捕捉目標(biāo)密度的某個(gè)特征（可能導(dǎo)致不可忽略的偏差）時(shí)，通常會(huì)表現(xiàn)為 發(fā)散的模擬過(guò)程 （diverging simulation），這種情況可以被檢測(cè)出來(lái)，并提示我們需要使用更小的步長(zhǎng)；

• 在軌跡終止前就達(dá)到了最大步數(shù)限制，說(shuō)明探索效率較低；

• 缺失信息的貝葉斯比例

  （Bayesian Fraction of Missing Information，BMFI）[9] 可以量化動(dòng)量重新采樣的效果與邊緣能量分布之間的匹配程度，可用于檢測(cè)預(yù)熱階段適應(yīng)不良或探索效率低下的情況。

Stan [21] 是目前最流行的通用貝葉斯推斷軟件。Stan 使用 C++ 編寫，具有獨(dú)立的命令行接口，同時(shí)也提供了成熟且廣泛使用的 Python 和 R 接口 （分別是 RStan [122] 和 PyStan [101]），以及輕量級(jí)封裝（如 CmdStanPy [123]、CmdStanR [42] 和 BridgeStan 2）。

此外，Stan 還支持大多數(shù)傳統(tǒng)用于數(shù)據(jù)分析的語(yǔ)言，包括：

• Matlab（Matlabstan）

• Julia（Stan.jl）

• Stata（StataStan）

• Mathematica（MathematicaStan）

• Scala（ScalaStan）

• 基于 HTTP 請(qǐng)求的接口（httpstan）

盡管 Stan 也實(shí)現(xiàn)了黑盒變分推斷 [69]、拉普拉斯近似和標(biāo)準(zhǔn)優(yōu)化方法，其核心的貝葉斯計(jì)算方法仍然是 NUTS （No-U-Turn Sampler），即 HMC 的一種變體。Stan 擁有豐富的數(shù)學(xué)庫(kù)，并集成了自動(dòng)微分（AD）[22]，還支持基于 OpenCL 的 GPU 加速與核融合（kernel fusion）[23, 24]。

Stan 的概率編程語(yǔ)言

Stan 的 PPL 是一種命令式語(yǔ)言 ，用戶用它來(lái)指定（對(duì)數(shù)）后驗(yàn)分布的計(jì)算過(guò)程。一個(gè) Stan 程序被劃分為多個(gè)代碼塊，其中最重要的是：

• data：輸入數(shù)據(jù)；
• parameters：待估計(jì)參數(shù)；
• model：模型定義。

參見附錄中提供的一個(gè)用 Stan 編寫的模型示例。

在 Stan 中，數(shù)據(jù)與參數(shù)的區(qū)別是在編譯時(shí)就確定的 ，因此如果要將某個(gè)變量從數(shù)據(jù)改為參數(shù)（或反過(guò)來(lái)），需要將其從 data塊移到 parameters 塊（或相反），并重新編譯程序。

關(guān)于 Stan 語(yǔ)言的重要研究工作包括 SlicStan [51, 52, 53]，它引入了多項(xiàng)改進(jìn)；還有將 Stan 轉(zhuǎn)換為 Pyro 的嘗試 [4]。

由于 Stan 基于 HMC 進(jìn)行計(jì)算，它所能擬合的模型類別是那些具有連續(xù)密度函數(shù) 的模型。一個(gè)重要的限制是：不能直接處理含有離散參數(shù)的模型 ，這些模型目前必須通過(guò)手動(dòng)邊緣化（marginalization）進(jìn)行處理。

這意味著 Stan 并不能完全替代 BUGS 所能實(shí)現(xiàn)的功能，HMC 也不會(huì)讓基于 Gibbs 抽樣的軟件變得過(guò)時(shí)。然而，實(shí)證研究表明，在適用的情況下，Stan 目前仍是通用貝葉斯推斷的首選工具 [7]。

值得注意的是，大多數(shù) Stan 用戶并不會(huì)直接編寫 Stan 語(yǔ)言代碼 。目前已有一些流行的數(shù)據(jù)分析包提供了簡(jiǎn)化版的公式或選項(xiàng)式建模語(yǔ)言，并使用 Stan 作為后臺(tái)進(jìn)行建模與計(jì)算：

• brms [16]：R 包，用于擬合層次模型（hierarchical models）；
• Prophet [117]：分別以 Python [119] 和 R [118] 實(shí)現(xiàn)，用于帶有趨勢(shì)、季節(jié)性和節(jié)假日效應(yīng)的非線性時(shí)間序列預(yù)測(cè)；
• rstanarm [6]:R包，提供R中lm、glm、aov等函數(shù)的貝葉斯版本。

總體而言，已有超過(guò) 140 個(gè)基于 Stan 的 R 包 ，它們?yōu)楦鞣N常見應(yīng)用領(lǐng)域中的典型模型提供了簡(jiǎn)單易用的接口。這些包的成功不僅歸功于 Stan 引擎本身，還得益于 R、Python 和 Julia 中日益豐富的實(shí)用工具的發(fā)展。

3.2 ADMB

AD Model Builder（ADMB） [38] 是第一個(gè)基于自動(dòng)微分（AD）的概率編程語(yǔ)言。它與 Stan 類似，數(shù)據(jù)、參數(shù)、似然函數(shù)和先驗(yàn)分布都是分開描述的，并且在編譯時(shí)就區(qū)分了數(shù)據(jù)與參數(shù)。

ADMB 更側(cè)重于基于優(yōu)化的推斷，但也實(shí)現(xiàn)了 Metropolis-Hastings（MH）、拉普拉斯近似，以及需手動(dòng)調(diào)參的 HMC。目前也有第三方實(shí)現(xiàn)的 NUTS 算法用于 ADMB [84]。

3.3 PyMC

PyMC [105] 是一個(gè)用于貝葉斯推斷的 Python 庫(kù)。它包括 HMC、SMC 和黑盒變分推斷方法。

PyMC 基于 PyTensor [^3]，這是一個(gè) Python 數(shù)學(xué)庫(kù)，是 Aesara 的一個(gè)分支，也是曾經(jīng)流行的 Theano 庫(kù)的延續(xù)。Theano 曾作為 PyMC3 的后端，直到當(dāng)前主要版本 4.0 發(fā)布并更名為 PyMC。

[^3]: 原文為 "PyTensor, 3"，推測(cè)為排版錯(cuò)誤，應(yīng)為 "PyTensor, a Python mathematics library..."

PyMC 中的計(jì)算圖可以被轉(zhuǎn)換為 C 代碼、Numba 或 JAX [15]（一種高性能的自動(dòng)微分庫(kù)，可將 Python/NumPy 代碼運(yùn)行在 CPU、GPU 和 TPU 上），這使得生成的代碼高度優(yōu)化。

PyMC 的語(yǔ)法與其他貝葉斯軟件中的概率編程語(yǔ)言相似。附錄中提供了一個(gè)用 PyMC 編寫的模型示例。

此外，還有一個(gè)基于 PyMC 構(gòu)建的包叫做 Bambi（BAyesian Model-Building Interface） [19]，它旨在簡(jiǎn)化層次廣義線性模型（GLMs）的使用體驗(yàn)。

4. 其他軟件 4.1 R-INLA

R-INLA [71, 130] 是一個(gè)流行的用于貝葉斯推斷的 R 包，專為**潛在高斯模型（latent Gaussian models）**設(shè)計(jì)。這類模型不需要使用概率編程語(yǔ)言（PPL），而是通過(guò)標(biāo)準(zhǔn)的 R 公式語(yǔ)法來(lái)指定模型，類似于 R 核心統(tǒng)計(jì)包 stats 中的 lm / glm 函數(shù)，并擴(kuò)展了用于平滑項(xiàng)和層次結(jié)構(gòu)（“隨機(jī)效應(yīng)”）的公式語(yǔ)法，例如在 R 包 mgcv [137] 中使用的語(yǔ)法。

當(dāng)超參數(shù)數(shù)量適中并滿足一些附加假設(shè)時(shí)，潛在高斯模型可以使用 集成嵌套拉普拉斯近似（INLA） [82, 103, 104] 這一高效的近似貝葉斯計(jì)算方法進(jìn)行快速計(jì)算。對(duì)于符合這些條件的模型來(lái)說(shuō)，R-INLA 是一種非常高效的 MCMC 替代方案，難以被取代。

R-INLA 的一個(gè)關(guān)鍵特性是支持使用隨機(jī)偏微分方程（SPDE） 方法建模連續(xù)空間模型 [3, 68, 72]。

最近，R-INLA 在模型表示方面進(jìn)行了改進(jìn)，并用變分貝葉斯修正層替換了原有的拉普拉斯近似，以提升其在數(shù)據(jù)規(guī)模、模型規(guī)模以及計(jì)算核心數(shù)量方面的可擴(kuò)展性 [44, 131, 132]。這一改進(jìn)還可以擴(kuò)展到對(duì)超參數(shù)的方差、偏度建模，以及對(duì)邊緣分布進(jìn)行修正。

4.2 通用 PPLs（Universal PPLs）

目前還沒有關(guān)于“通用 PPL”的統(tǒng)一定義。但普遍接受的是：一個(gè)通用 PPL 程序可以在任何地方包含概率操作，例如，甚至無(wú)法靜態(tài)確定隨機(jī)變量的數(shù)量。我們將這種能力稱為 [100] 中提到的“反轉(zhuǎn)模擬器”（inverting simulators）的概念：用戶編寫一個(gè)隨機(jī)模擬程序，而 PPL 框架可以根據(jù)觀測(cè)數(shù)據(jù)推斷出該模擬程序的性質(zhì)。

從這個(gè)意義上講，BUGS、Stan 和前面提到的其他語(yǔ)言并不是通用 PPLs ；它們更像是為特定模型類別優(yōu)化的貝葉斯推斷系統(tǒng)，使得這些類別的模型可以自動(dòng)進(jìn)行推斷。而設(shè)計(jì)一個(gè)通用 PPL 的核心目標(biāo)是構(gòu)建一個(gè)通用編程語(yǔ)言 ，并配有一個(gè)能夠處理所有可能算法的推理框架。

理論上，通用 PPL 可以涵蓋整個(gè)貝葉斯推斷過(guò)程，而且編寫一個(gè)隨機(jī)模擬器通常比設(shè)計(jì)合適的統(tǒng)計(jì)推斷方法更容易。然而，是否可以實(shí)現(xiàn)對(duì)如此廣泛類別的算法進(jìn)行高效且自動(dòng)化的推理 ，目前仍不明確。

從貝葉斯統(tǒng)計(jì)實(shí)踐者的角度來(lái)看，通用 PPLs 目前仍然是研究對(duì)象，而非廣泛實(shí)用的工具。不過(guò)，已經(jīng)出現(xiàn)了許多令人振奮的發(fā)展。在本節(jié)剩余部分，我們將介紹一些較為流行或較新的通用 PPLs。

其他相關(guān)的早期通用 PPL 語(yǔ)言包括：

• Church [50]

• Venture [79]

• Anglican [128]

• 基于 Julia 的 Gen [28]

• Turing.jl [45]（及其更近期的前端 DynamicPPL [116]）

• 嵌入 Python 的 Edward/Edward2 [129]

Bean Machine 是一個(gè)基于貝葉斯的軟件，也是一種聲明式的通用 PPL ，嵌入在 Python 中，后端使用 PyTorch。

本質(zhì)上，Bean Machine 支持對(duì)具有不同變量數(shù)量的貝葉斯網(wǎng)絡(luò)進(jìn)行分布建模。盡管命令式語(yǔ)言（如當(dāng)前最流行的 Stan）正變得越來(lái)越常見，作者仍然主張采用聲明式 PPL。他們認(rèn)為，變量之間的（可能是動(dòng)態(tài)的）依賴關(guān)系更容易從聲明式模型描述中恢復(fù)出來(lái)，并且推理可以更靈活地適應(yīng)不同的變量塊，包括在高維情況下通常不可行的二階方法。

Bean Machine 實(shí)現(xiàn)了多種單點(diǎn)采樣器、NUTS、牛頓蒙特卡洛（Newtonian Monte Carlo）和黑盒變分推斷（black-box VI）。它還支持變量分塊（blocking）和自定義提議器（custom proposers）。

附錄中提供了一個(gè)用 Bean Machine 編寫的模型示例。

4.2.2 Birch [88]

Birch 是一個(gè)通用 PPL，它將代碼轉(zhuǎn)換為 C++，并支持 GPU 加速。用戶以生成式方式實(shí)現(xiàn)模型的聯(lián)合分布，偏好泛型與面向?qū)ο缶幊谭妒健?/p>

Birch 的推理方法基于序貫蒙特卡洛（SMC） ，并結(jié)合了基于梯度的核函數(shù)（kernel-based）方法。

Birch 的一大特色是支持自動(dòng)邊緣化 （automatic marginalization）和自動(dòng)條件設(shè)定 （automatic conditioning）。與自動(dòng)微分類似，Birch 能識(shí)別已知形式（如共軛先驗(yàn)和離散枚舉），從而在可能的情況下對(duì)隨機(jī)變量進(jìn)行邊緣化，并在必要時(shí)在后續(xù)模擬中進(jìn)行條件設(shè)定。

這些功能的實(shí)現(xiàn)基于一種被稱為延遲抽樣 （delayed sampling）[87] 的啟發(fā)式方法，在程序執(zhí)行過(guò)程中盡可能推遲隨機(jī)變量的模擬，從而揭示出可以利用的機(jī)會(huì)。

結(jié)果是，推理方法得到了諸如 Rao-Blackwellization [87] 和 變量消除 （variable elimination）[136] 等增強(qiáng)功能。

Birch 已經(jīng)成功應(yīng)用于一些未知隨機(jī)變量數(shù)量的問(wèn)題中，例如：

• 多目標(biāo)跟蹤 [88]（物體數(shù)量未知）；

• 統(tǒng)計(jì)系統(tǒng)發(fā)育學(xué) [102]（系統(tǒng)發(fā)育樹中的滅絕分支數(shù)量未知）。

附錄中提供了一個(gè)用 Birch 編寫的模型示例。

4.2.3 Pyro [12]

Pyro 是一個(gè)基于 Python 的概率編程語(yǔ)言（PPL），構(gòu)建在 PyTorch [93] 后端之上。Pyro 的主要計(jì)算方法是隨機(jī)變分推斷（stochastic variational inference），因此它主要面向可擴(kuò)展的概率機(jī)器學(xué)習(xí)任務(wù)。

NumPyro [94] 是 Pyro 的一個(gè)基于 NumPy 的后端，使用 JAX 進(jìn)行自動(dòng)微分，并支持編譯到 CPU/GPU 上運(yùn)行。

4.2.4 Blang [14]

Blang 是一個(gè)開源軟件包，用于對(duì)任意空間上的后驗(yàn)分布進(jìn)行近似，即支持不僅包含整數(shù)和實(shí)數(shù)變量，還包含用戶自定義數(shù)據(jù)類型的貝葉斯模型，例如系統(tǒng)發(fā)育樹、隨機(jī)圖和序列比對(duì)等。

Blang 項(xiàng)目包括一個(gè)標(biāo)準(zhǔn)庫(kù)，其中包含了用 Blang 語(yǔ)言編寫的一些常見數(shù)據(jù)類型和分布，以及用于創(chuàng)建新數(shù)據(jù)類型及其相關(guān)分布的擴(kuò)展接口。

用戶可以發(fā)布包含新數(shù)據(jù)類型和分布的版本化 Blang 包，并導(dǎo)入他人貢獻(xiàn)的包及其傳遞依賴項(xiàng)。

Blang 語(yǔ)言的作用域規(guī)則被用于自動(dòng)檢測(cè)稀疏模式，并構(gòu)建一種稱為因子圖 （factor graph）的圖形模型。基于這個(gè)因子圖，后驗(yàn)分布通過(guò)一種自適應(yīng)不可逆并行退火算法 （adaptive non-reversible parallel tempering algorithm）[114] 來(lái)近似，默認(rèn)情況下在用戶的 CPU 核心上并行執(zhí)行，也可以通過(guò)與 Pigeons 分布式并行退火包集成，在 MPI（消息傳遞接口）環(huán)境中分布式執(zhí)行。

附錄中提供了一個(gè)用 Blang 編寫的模型示例。

4.3 無(wú)需似然的貝葉斯推斷（Likelihood-free Bayesian Inference）

無(wú)需似然的推斷方法 （Likelihood-free inference, LFI），如：

• 近似貝葉斯計(jì)算

  （Approximate Bayesian Computation, ABC）[108]

• 貝葉斯合成似然

  （Bayesian Synthetic Likelihood, BSL）[99]

• 基于機(jī)器學(xué)習(xí)的后驗(yàn)近似

• 替代似然方法 [56, 25]

這些方法適用于無(wú)法或難以計(jì)算似然函數(shù)，但存在生成式模擬器模型的情況。這類方法廣泛應(yīng)用于天體統(tǒng)計(jì)學(xué)、遺傳學(xué)、生態(tài)學(xué)、系統(tǒng)生物學(xué)和人類認(rèn)知建模等領(lǐng)域。

ELFI（Engine for Likelihood-Free Inference） [73] 是一個(gè)用于 LFI 的 Python 軟件包，涵蓋了上述所有主要方法（ABC、BSL 和基于 ML 的方法）。ELFI 具有模塊化設(shè)計(jì)，由一個(gè)基于 DAG 的建模 API 和一個(gè)獨(dú)立的推理 API 構(gòu)成，允許用戶靈活選擇各種算法來(lái)從近似后驗(yàn)分布中抽樣。

抽樣可以通過(guò)自適應(yīng)重要性抽樣（Adaptive Importance Sampling）或 MCMC/HMC 實(shí)現(xiàn)，并可以選擇是否使用替代模型來(lái)近似似然函數(shù)。替代模型使用高斯過(guò)程（GPs）和主動(dòng)學(xué)習(xí)（貝葉斯優(yōu)化）來(lái)模擬目標(biāo)函數(shù)。主動(dòng)學(xué)習(xí)方法已被證明可以將無(wú)需似然的推斷加速幾個(gè)數(shù)量級(jí)。

其他通用的 ABC 軟件包包括：

• Python：pyABC [106]、ABCpy [34]

• R：abc [26]、EasyABC [64]、abctools [92]、ABCreg [127]

基于神經(jīng)網(wǎng)絡(luò)的替代模型可通過(guò) Python 包 sbi [125] 和 R 包 [2] 提供的 BSL 工具箱實(shí)現(xiàn)。

關(guān)于 ABC 軟件的更詳細(xì)綜述見 Nunes 和 Prangle [92] 以及 Kousathanas 等人 [67]。

4.4 以計(jì)算為核心的軟件

• Blackjax ：一個(gè)為 JAX 設(shè)計(jì)的 Python MCMC 方法庫(kù)。它可在 CPU 和 GPU 上運(yùn)行，具有魯棒性和高效性，并能輕松集成到支持 JAX 兼容密度函數(shù)的 PPL 中（如 TFP、Oryx、NumPyro、Aesara、PyTensor/PyMC）。

• emcee [36, 37]：一個(gè) Python 實(shí)現(xiàn)的仿射不變 MCMC 集合采樣方法（Affine Invariant MCMC Ensemble Sampler）[49]。這種無(wú)需導(dǎo)數(shù)的方法適用于低維問(wèn)題中的黑盒似然函數(shù)，這在天體物理學(xué)中非常常見。

• dynesty [109]：另一個(gè)在天體物理學(xué)中流行的 Python 包，實(shí)現(xiàn)了動(dòng)態(tài)嵌套采樣 （dynamic nested sampling）[58]。

• Mamba.jl [^6]：一個(gè)面向希望使用和開發(fā) MCMC 方法的用戶的 Julia 包。它實(shí)現(xiàn)了多種 MCMC 方法（如 HMC、NUTS、Metropolis-within-Gibbs 等）以及 MCMC 診斷工具。

• DynamicHMC.jl ：另一個(gè)流行的 Julia 包，實(shí)現(xiàn)了最先進(jìn)的貝葉斯計(jì)算方法。

[^6]: GitHub 地址：https://github.com/brian-j-smith/Mamba.jl

4.5 其他通用軟件

• Infer.NET [83]：一個(gè)用 C# 編寫的機(jī)器學(xué)習(xí)庫(kù)，適用于 .NET 框架。它支持對(duì)貝葉斯網(wǎng)絡(luò)和馬爾可夫隨機(jī)場(chǎng)進(jìn)行自動(dòng)近似推斷。貝葉斯計(jì)算主要限于消息傳遞方法。

• TensorFlow Probability (TFP) [31]：一個(gè)構(gòu)建在 TensorFlow [35] 之上的 Python 庫(kù)。

• greta [48]：一個(gè)語(yǔ)法非常簡(jiǎn)潔的概率編程語(yǔ)言，嵌入在 R 中，基于 TensorFlow 和 TFP。雖然其提供的貝葉斯計(jì)算方法有限，但它具有良好的可擴(kuò)展性。附錄中提供了一個(gè)用 greta 編寫的模型示例。

• Oryx ^8 ：一個(gè)構(gòu)建在 JAX 之上的 PPL。

• 《統(tǒng)計(jì)軟件雜志》 （Journal of Statistical Software）最近還出版了一期關(guān)于貝葉斯軟件的專刊 [18]，其中包括本文介紹的一些軟件以及其他專門用途的軟件。

CODA [96] 是一個(gè)至今仍廣泛使用的 R 包，用于對(duì) MCMC 輸出進(jìn)行事后診斷和分析。

ArviZ [70] 是一個(gè) Python 包，提供 MCMC 診斷、模型評(píng)估和模型驗(yàn)證工具。ArviZ 不依賴于后端計(jì)算框架，是目前 Python 中最受歡迎的此類工具。

bridgesampling （R 包）[55] 可以估計(jì)邊緣似然、貝葉斯因子、后驗(yàn)?zāi)Ｐ透怕室约皻w一化常數(shù)。

posterior （R 包）[17] 提供了對(duì) MCMC 樣本的子集選取、綁定、變換以及格式轉(zhuǎn)換功能，并包含了一些輕量級(jí)實(shí)現(xiàn)的最先進(jìn)后驗(yàn)推斷診斷工具。

BayesFactor [85] 是一個(gè)用于計(jì)算貝葉斯因子的 R 包，適用于列聯(lián)表、單樣本與雙樣本設(shè)計(jì)、單因素設(shè)計(jì)、一般方差分析（ANOVA）設(shè)計(jì)以及線性回歸模型。

shinystan （R 包）[43] 提供了一個(gè)圖形用戶界面，用于交互式地進(jìn)行馬爾可夫鏈蒙特卡洛（MCMC）診斷及其他用于分析后驗(yàn)抽樣結(jié)果的工具。shinystan 對(duì) MCMC 樣本的生成方式不敏感，但對(duì)使用 RStan 擬合的模型具有一些額外的功能支持。

loo （R 包）[133] 為使用 MCMC 方法擬合的貝葉斯模型提供高效的近似留一交叉驗(yàn)證（leave-one-out cross-validation）方法。

bayestestR （R 包）[78] 提供了一套在貝葉斯統(tǒng)計(jì)框架下處理不確定性和效應(yīng)的工具。它對(duì)后驗(yàn)樣本的生成軟件不敏感，包括最大后驗(yàn)估計(jì)（MAP）、離散度度量、ROPE 分析（Region of Practical Equivalence）以及貝葉斯因子等。

bayesplot （R 包）[40] 提供了用于繪制貝葉斯模型的函數(shù)，包括后驗(yàn)抽樣結(jié)果、可視化 MCMC 診斷圖以及圖形化的預(yù)測(cè)檢查。

projpred （R 包）[95] 針對(duì)使用 MCMC 方法擬合的貝葉斯廣義線性模型和加性多層次模型，執(zhí)行投影預(yù)測(cè)變量選擇（projection predictive variable selection）。

priorsense （R 包）[66] 支持對(duì)使用 MCMC 方法擬合的貝葉斯模型進(jìn)行先驗(yàn)分布和似然函數(shù)的敏感性分析。

MCMCpack [80] 是一個(gè)實(shí)現(xiàn)基于 MCMC 的多種統(tǒng)計(jì)方法的 R 包。

用于貝葉斯網(wǎng)絡(luò)結(jié)構(gòu)學(xué)習(xí)和參數(shù)估計(jì)的工具包括：

• bnlearn [107]
• Matlab 的 Bayes Net 工具箱 [86]
• HUGIN [77]
• VIBES [13]
• MSBNx [65]

• 商業(yè)工具如 GEnIe/SMILE [32] 和 Netica

貝葉斯推斷軟件領(lǐng)域從未像今天這樣活躍和多樣化。各個(gè)方向都在不斷發(fā)展，不僅提供了更好的工具，使得典型建模任務(wù)更加易于使用、更具魯棒性或更高效，同時(shí)也在不斷拓展貝葉斯方法能夠?qū)崿F(xiàn)的邊界。

與編程語(yǔ)言類似，人們可能會(huì)根據(jù)任務(wù)選擇不同的語(yǔ)言：比如用 Python 進(jìn)行通用編程，用 R 進(jìn)行數(shù)據(jù)清洗和可視化，或者使用新興的 Julia 進(jìn)行高性能數(shù)據(jù)分析。對(duì)于貝葉斯推斷軟件來(lái)說(shuō)，也沒有一種“萬(wàn)能”的解決方案 。

• Stan 通常是貝葉斯建模與推斷的首選工具；
• Pyro 或 TFP（TensorFlow Probability） 更適用于貝葉斯機(jī)器學(xué)習(xí)；

• 而還有許多其他工具專門用于更特定的任務(wù)。

這種多樣性是可以理解的，因?yàn)橄拗乒ぞ叩墓δ芸梢院?jiǎn)化其設(shè)計(jì)，并允許更高效的計(jì)算。雖然在通用概率編程語(yǔ)言（PPL）和基礎(chǔ)貝葉斯計(jì)算方面已經(jīng)取得了一些令人鼓舞的進(jìn)展，但目前尚不清楚是否能夠在表達(dá)能力與效率之間找到新的平衡點(diǎn) ，從而催生第三代工具。

因此，用戶要么必須接受所選工具的局限性，要么學(xué)會(huì)使用多種工具和語(yǔ)言。一個(gè)自然的解決方案是實(shí)現(xiàn)語(yǔ)言之間的自動(dòng)轉(zhuǎn)換 ，或從統(tǒng)計(jì)表示法自動(dòng)生成代碼，如附錄 A 所示。然而這是一個(gè)極具挑戰(zhàn)性的問(wèn)題，因?yàn)椴煌Z(yǔ)言在表達(dá)能力上存在差異，即使存在自動(dòng)轉(zhuǎn)換機(jī)制，也可能生成效率低下的代碼。無(wú)論如何，這個(gè)領(lǐng)域似乎缺乏足夠的激勵(lì)機(jī)制來(lái)推動(dòng)這一方向的發(fā)展。

一個(gè)新的 PPL 往往是通過(guò)帶有代碼和數(shù)據(jù)的模型示例來(lái)學(xué)習(xí)的，或者是通過(guò)從我們已熟悉的語(yǔ)言進(jìn)行翻譯來(lái)學(xué)習(xí)的。因此，像 BUGS 和 Stan 這樣流行的 PPL 擁有詳盡的文檔就不足為奇了，包括用戶手冊(cè)、案例研究，以及在 Stan 的情況下，還提供了從 BUGS 到 Stan 的翻譯示例。

流行的 PPL 通常具備以下特征：

• 可在所有主流平臺(tái)上使用，并提供多種主要編程語(yǔ)言的接口（通常為獨(dú)立程序并配有接口）；

• 具有開放的治理結(jié)構(gòu)；

• 擁有一個(gè)活躍的社區(qū)，促進(jìn)用戶與開發(fā)者之間的交流，這通常從項(xiàng)目誕生之初就已經(jīng)建立。

隨著統(tǒng)計(jì)分析標(biāo)準(zhǔn)的不斷提高，對(duì)統(tǒng)計(jì)工作流程的支持 也變得越來(lái)越重要。某種程度上，一些優(yōu)秀的獨(dú)立工具已經(jīng)對(duì)此進(jìn)行了覆蓋。然而，工作流程中的某些部分較難封裝，因?yàn)樗鼈円蕾囉诘讓诱Z(yǔ)言或計(jì)算機(jī)制。

另一個(gè)問(wèn)題是，當(dāng)使用黑盒組件時(shí)，如果分析結(jié)果不如預(yù)期，從軟件角度來(lái)看有時(shí)很難定位問(wèn)題所在。

現(xiàn)代貝葉斯推斷軟件已經(jīng)普遍認(rèn)識(shí)到可擴(kuò)展性 （scalability）的實(shí)際重要性，幾乎所有流行的語(yǔ)言都至少通過(guò)第三方或原生庫(kù)提供了一定程度的支持。

• 針對(duì) 數(shù)據(jù)規(guī)模的可擴(kuò)展性 ，當(dāng)前主要依賴于 GPU 等大規(guī)模并行設(shè)備上的優(yōu)化矩陣代數(shù)計(jì)算。我們預(yù)計(jì)，隨著硬件和計(jì)算庫(kù)的發(fā)展，這方面的能力還將進(jìn)一步提升。

• 而針對(duì) 模型規(guī)模的可擴(kuò)展性 則更多地取決于所使用的貝葉斯計(jì)算方法。目前， HMC 是通用全貝葉斯計(jì)算的最先進(jìn)方法 。

盡管 HMC 存在一些局限性（與其他通用貝葉斯計(jì)算方法類似），但這些問(wèn)題通常可以通過(guò)精心的重新參數(shù)化 ，或?qū)⒉煌挠?jì)算任務(wù)分配給不同的參數(shù)塊來(lái)緩解。然而，目前還沒有自動(dòng)化的方式來(lái)實(shí)現(xiàn)這一點(diǎn)。

未來(lái)，通用軟件（如 Stan、JAGS、MultiBUGS）將如何被取代，以及通用 PPL、近似貝葉斯計(jì)算（ABC）和序貫蒙特卡洛（SMC） 將扮演什么角色，仍有待觀察。

附錄 A：建模語(yǔ)言示例

在本附錄中，我們將通過(guò)一個(gè)貝葉斯線性回歸 的示例來(lái)展示幾種建模語(yǔ)言的具體寫法。

原文鏈接： https://sites.stat.columbia.edu/gelman/research/published/Bayesian_software_review-8.pdf