校準貝葉斯推理

Calibrating Bayesian Inference

校準貝葉斯推理

摘要
盡管貝葉斯統計因其直觀的不確定性量化和靈活的決策能力而在心理學研究中廣受歡迎，但其在有限樣本中的表現可能不可靠。本文揭示了一個關鍵弱點：當研究者所選的先驗分布與真實參數生成過程不匹配時，貝葉斯推斷在長期運行中可能產生誤導。鑒于真實過程在實踐中通常未知，我們提出一種更安全的替代方案：對貝葉斯可信區域進行校準，以實現頻率學派的有效性。這一標準更強，能保證無論底層參數生成機制如何，貝葉斯推斷均有效。為在實踐中解決校準問題，我們提出了一種新穎的隨機逼近算法。我們開展并報告了一項蒙特卡洛實驗，結果表明，在某些參數生成情形下，未經校準的貝葉斯推斷可能過于寬松（liberal），而我們的校準方案始終能保持有效性。

關鍵詞：貝葉斯推斷，頻率學派推斷，統計有效性，可信區域，隨機逼近，黎曼優化

引言

近幾十年來，使用貝葉斯方法的心理學出版物顯著增加（例如，Kruschke, 2021；van de Schoot 等, 2021；van de Schoot, Winter, Ryan, Zondervan-Zwijnenburg & Depaoli, 2017；Volpe 等, 2025）。貝葉斯推斷通過后驗概率度量提供直觀的不確定性量化。借助從后驗分布中抽取的（近似）隨機樣本，可方便地以無需解析推導的方式對模型參數進行推斷、預測未來數據并評估模型擬合（Gelman, Carlin, Stern & Rubin, 2013）。

現成的貝葉斯分析軟件不僅包括通用的馬爾可夫鏈蒙特卡洛（MCMC）抽樣器，如 JAGS（Plummer, 2017）和 Stan（Stan Development Team, 2024），還包括專為特定建模框架設計的程序，如 Mplus（Muthén & Muthén, 1998–2024）和 blavaan（Merkle & Rosseel, 2018）。

貝葉斯推斷的哲學與統計基礎在于在不確定性下做出一致的決策（DeGroot, 1970；Savage, 1954）：用恰當的概率測度編碼先驗知識，并在觀測新數據后通過貝葉斯公式更新知識。然而，貝葉斯方法在心理學中的應用大多是工具性的，而非哲學性的。將先驗明確聯系到研究者信念的正當性說明仍然罕見，而依賴默認或共軛先驗的做法卻十分普遍。通常有三種理由被用來為此辯護。

第一，某些默認先驗——尤其是“弱信息”和“客觀”先驗（例如 Berger, Bernardo & Sun, 2024, 2015；Datta & Mukerjee, 2004；Gelman, Jakulin, Pittau & Su, 2008）——在現有文獻中已被證明具有良好的理論性質或強健的實證表現。第二，在適當的正則條件下，隨著樣本量增大，先驗的影響會減弱，所得的后驗推斷在大樣本下通常與頻率學派結果相似（例如 Bernstein-von Mises 定理；van der Vaart, 1998，第10章）。第三，通常建議進行敏感性分析，以確保統計結論對不同先驗選擇的穩健性（例如 Depaoli, 2022；Depaoli, Winter & Visser, 2020；Van Erp, Mulder & Oberski, 2018）。

然而，當應用于現實世界的心理學研究時，上述辯護的合理性常常值得懷疑。首先，“客觀性”和“無信息性”在定義默認先驗時并未建立在單一、統一的框架之上（例如 Kass & Wasserman, 1996）。此外，哪種默認先驗表現最佳往往取決于具體模型（Yang & Berger, 1998）。因此，尋找一個在所有情況下都表現優異的先驗很可能是一個難以實現的目標。

第二，由于研究焦點或實際考量，實質性研究者可能不得不處理小樣本。例如，由多重社會身份定義的交叉性亞群體通常過于狹窄，難以積累足夠數據（例如 Cole, 2009）。基于大樣本理論的統計程序在小樣本應用中可能數值不穩定，甚至產生誤導性推斷（例如 van de Schoot & Mio?evi?, 2020）。

第三，先驗敏感性分析可能無法得出明確結論。幾乎總能找到一種病態的先驗分布（例如，其質量集中在遠離原始貝葉斯解的區域），從而推翻原有結論。此外，貝葉斯計算可能計算成本過高，難以重復大量次數。因此，先驗敏感性分析通常局限于一組有限且任意選擇的先驗，對先驗設定的診斷價值甚微。

為應對現實中“實用主義貝葉斯”（pragmatic Bayes）的普遍性，貝葉斯方法的方法論研究越來越關注在數據和/或參數重復抽樣下的表現。事實上，任何在長期運行中系統性地產生錯誤結論的推斷程序都應被摒棄。為此，過去幾十年開展了大量蒙特卡洛（MC）實驗，其中貝葉斯推斷程序（即假設檢驗和區間估計）在各種數據與參數生成機制及設計因素（如樣本量、協變量數量等）下被評估（例如 Finch & French, 2019；McNeish, 2016, 2017a, 2017b；Preacher & MacCallum, 2002；Smid, McNeish, Mio?evi? & van de Schoot, 2020）。然而，MC 研究的一個主要局限在于其結論依賴于特定模型和設計，難以推廣到未明確測試的情境之外。因此，目前使用貝葉斯分析的心理學研究所得結論的可信度尚未完全確立。

本文有兩個主要目標：一是教學性的，二是方法論的。借鑒統計決策理論（Berger, 1985）和推斷模型（Inferential Models, IMs；C. Liu & Martin, 2024；Martin & Liu, 2015）的相關成果，我們強調貝葉斯有效性（Bayesian validity）這一核心概念：即關于參數的不合理陳述不太可能頻繁地具有較大的后驗概率。這一原則為貝葉斯程序的長期頻率表現提供了正當性（例如 Martin, 2022a, 2022b, 2022c；精確的定義見“貝葉斯推斷與統計有效性”一節）。在缺乏先驗知識、僅將貝葉斯方法作為工具使用的場景下，我們指出應追求更強的頻率學派有效性（frequentist validity），因為它必然蘊含對任意先驗設定下的貝葉斯有效性。

在方法論層面，我們提出一種計算程序，用于校準基于后驗的推斷以確保頻率學派有效性。該方法結合了無梯度隨機逼近（gradient-free stochastic approximation, SA）與流形優化（manifold optimization）（例如 Absil, Mahony & Sepulchre, 2008；Spall, 1992）。我們在一項蒙特卡洛實驗中將該方法應用于線性-正態因子分析模型（例如 J?reskog, 1969），結果表明：若未經適當校準，貝葉斯推斷可能是無效的。

本文其余部分結構如下：我們首先回顧統計決策理論、推斷模型（IMs）和統計有效性的理論基礎。具體而言，我們闡明一個關鍵事實：頻率學派有效性可確保對任意先驗選擇下的貝葉斯有效性。接著，我們介紹一種實用的計算算法，用于校準貝葉斯推斷以實現頻率學派有效性。我們開展了一項概念驗證型蒙特卡洛（MC）實驗，將基于隨機逼近（SA）的校準推斷與基于MCMC抽樣的原始推斷進行對比。最后，文章通過討論本研究的啟示、局限性及未來研究方向予以總結。

貝葉斯推理與統計有效性

貝葉斯模型和嵌套置信區域

后驗可能性

嵌套置信區域不僅僅是模型參數的一組估計量。它們在貝葉斯推理中的基礎作用可以通過可能性理論正式建立（Dubois, 2006; Dubois & Prade, 1988; Zadeh, 1978）。設

后驗可能性輪廓可以通過在所有 α 級上拼接嵌套置信區域族來直觀地描繪；圖1中可以找到圖形說明。更深入的可能性理論及其在統計推理中的應用可以在例如Denoeux和Li (2018)，C. Liu和Martin (2024)，以及Martin (2025b)中找到。

后驗可能性是貝葉斯框架內推斷的基石。可能性等值線函數 ? 為所有嵌套可信區域提供了一個簡潔的總結：對于任意給定的 α ∈ [0,1]，可通過取該等值線的上 α-水平集直接獲得一個集合估計量。此外，在貝葉斯檢驗中，零假設的可信度由其后驗可能性保守地量化。例如，若簡單假設 H = {θ?} 的可能性 Π??{H} = ?(θ?) 小于預設的（較小的）α 水平，則該假設被拒絕。

統計有效性

為了使任何統計程序在實踐中可靠，最好能確立該程序在各種情境下均能產生可靠的結果。特別是在對模型參數進行推斷時，重要的是不應在長期運行中反復將低可能性值賦予頻繁發生的事件。越來越多的貝葉斯和頻率學派統計學家，盡管對概率代表什么持有不同觀點，但仍支持評估推斷程序在重復抽樣下表現的重要性（例如 Grünwald, 2018；Martin, 2022a, 2022b, 2022c）。接下來，我們將回顧貝葉斯有效性和頻率學派有效性的概念。我們強調，基于后驗可能性測度的貝葉斯推斷，在參數與數據模型均被正確設定的假設下，滿足貝葉斯有效性。同時，具有頻率學派有效性的程序，在數據模型被正確設定的前提下，對所有先驗而言自動滿足貝葉斯意義下的有效性。

貝葉斯有效性

令 h : Y × Q → ? 為任意 P_{Y,Θ}-可積函數。根據富比尼定理（Fubini’s）

定理（Billingsley, 2012，定理 18.3）指出，我們可以將函數 h 關于數據和參數的聯合期望寫成如下迭代期望的形式：

校準置信區域

在本節中，我們提出了一種通用的計算策略，通過閾值化觀測檢驗統計量來校準后驗可能性輪廓并確保頻率論有效性（8）。作為輸入，我們定義了一個通過閾值化觀測檢驗統計量來定義的嵌套置信區域族。然后，我們使用無梯度SA算法基于檢驗統計量的生存函數來校準置信區域。值得注意的是，我們解決的優化程序本質上等同于為IM可能性輪廓（9）找到一個“變分近似”（Cella & Martin, 2024; Martin, 2025a）。我們提案的一個獨特貢獻是同時擾動SA（Spall, 1992, 2000, 2009）和流形優化的新組合，這增強了校準算法的可擴展性，并促進了其在現實規模模型中的應用。

檢驗統計量和嵌套置信區域

具體在當前工作中，我們從兩種類型的檢驗統計量構建嵌套置信區域：Wald統計量和后驗密度比（PDR）統計量，分別生成橢圓形和HPD置信區域。

Wald 統計量與橢圓區域

與Wald統計量相關的嵌套置信區域族可以表示為：

后驗密度比統計量和最高后驗密度區域

或者，我們可以根據PDR統計量定義置信區域：

（13）是標準大樣本理論中最大似然估計的似然比統計量的推廣。還要注意，PDR統計量（13）是對Martin（2022b）在更一般的部分先驗背景下定義的“相對合理性排序”的對數變換。類似于（12），可以通過收集具有足夠低PDR統計量值的參數值來構建置信區域：

優化問題中，其可行區域形成歐幾里得空間的可微分子流形，可以通過黎曼梯度算法有效解決（Absil et al. 2008；另見Y. Liu 2020, 2021，關于黎曼梯度算法及其在心理測量問題中的應用的可訪問介紹）。特別地，黎曼梯度上升算法將歐幾里得空間中的常規梯度上升算法推廣，具有兩個顯著區別。首先，最陡上升方向只能在當前迭代處局部定義，通過黎曼梯度獲得，這是通過將目標函數的周圍梯度投影到子流形的切空間來實現的。對于我們的問題（15），

在應用標準黎曼梯度算法求解[15]時，另一個挑戰來自于對p值函數的精確梯度?θπy(θ)的評估。這是因為p值函數是一個積分，其定義域通過檢驗統計量依賴于θ。盡管有時可以通過廣義萊布尼茨法則（也稱為雷諾茲輸運定理；參見例如Flanders, 1973; Reddiger & Poirier, 2023）獲得解析導數，但一般而言，計算這些導數是困難的，因為檢驗統計量本身的計算可能涉及優化。一種簡單的嘗試是用有限差分（FD）估計來近似梯度的第r個元素（r = 1, ..., q），記為?θrπy(θ)：

對于某個小的擾動 c > 0，其中 er 是一個在第 r 個元素處為 1、其余位置為 0 的單位向量。由于方括號內的項在(20)中是其期望值的無偏估計量，我們或許可以期望通過黎曼隨機逼近（SA）算法（例如，Bonnabel, 2013; Tripuraneni, Flammarion, Bach, & Jordan, 2018）找到(15)的一個近似解。然而，這種近似存在兩個問題：(a) 如果 c 固定，有限差分近似的偏差不會消失；(b) 隨著參數數量增加，需要對檢驗統計量進行更多次評估，從而導致“維度災難”。

為解決這些問題，我們提出一種SPRSA算法，它結合了針對所有 q 個參數的同時擾動有限差分法（解決了問題(a)），以及沿黎曼梯度迭代 k = 0, 1, ... 進行衰減的有限差分步長序列 {ck}（解決了問題(b)）。設 Y = g(U, θ) 為一個數據生成算法，其中隨機分量 U ~ PU，且分布 PU 完全已知。在第 k 次迭代時，?θπy(θ(k)) 的同時擾動有限差分估計量定義為

使用SPRSA算法實現實際效率，需要仔細地、逐案調整多個方面：學習率序列{ak}、有限差分率序列{ck}以及總迭代次數K。在我們的數值研究中，SPRSA算法的調參細節將在“模擬研究”部分的“抽樣與調參細節”中提供。

蒙特卡洛實驗

數據生成

協方差結構模型已被廣泛用于解釋觀測響應變量之間的相關模式（Bollen, 1989; J?reskog, 1970）。令 Z = (Z?, ..., Z?)? ∈ ???? 為獨立同分布（i.i.d.）的樣本數據，其中每個 m 維響應向量 Z? ~ N(μ, Σ(θ))，其均值向量 μ ∈ ??，協方差矩陣 Σ(θ) ∈ ????? 由參數 θ ∈ ?q 參數化。在我們的數值示例中，我們考慮一個一維公共因子模型（J?reskog, 1969）：

在(24)中，λ = (λ?, ..., λ?)? ∈ ?? 是因子載荷向量，ψ ≥ 0 表示公共因子方差，u = (u?, ..., u?)? ∈ ??? 收集了唯一因子方差參數。為了識別該模型，要求公共因子 η? ∈ ? 與唯一因子 ε? ∈ ?? 相互獨立，并且第一個因子載荷 λ? 固定為1。以下協方差矩陣由公共因子模型推導得出：

為了避免對參數空間施加邊界限制，我們對所有方差參數進行了對數變換：

令 ζ = log ψ / 2 表示公共因子方差的對數標準差（SD），并且

ω? = log u? / 2，j = 1, ..., m，表示第 j 個唯一因子方差的對數標準差。我們將這些無界參數收集在參數向量 θ = (ζ, λ?, ..., λ?, ω?, ..., ω?)? 中，其維度為 q = 2m。由于 Z?, i = 1, ..., n 的獨立同分布正態性，樣本交叉乘積矩陣 Y = Σ???? (Z? ? Z?)(Z? ? Z?)?（其中 Z? 表示樣本均值向量）服從 Wish(Σ(θ), n?1) 分布，即尺度矩陣為 Σ(θ)、自由度為 n?1 的威沙特分布。由于交叉乘積矩陣 Y 是協方差結構 Σ(θ) 的充分統計量，我們在整個模擬研究中將 Y 視為數據。Y 的直接數據生成算法為 Y = g(U, θ) = Σ(θ)1/2 U Σ(θ)1/2，其中 U ~ Wish(I???, n?1)。

模擬條件由兩個交叉因素決定：響應變量的數量（m = 5 和 15），以及三種參數生成情景（情景1–3）。樣本量固定為 n = 100。在情景1中，響應變量的共性在各次重復中隨機從 U[.2, .8] 中抽取，涵蓋從低到高的共性水平（MacCallum, Widaman, Zhang, & Hong, 1999）。公共因子方差被設定為第一個響應變量的共性值。唯一方差被確定為使得所有響應變量具有單位方差。由此產生的 Θ 的分布作為 P*Θ。在情景2中，所有響應變量具有固定的低共性（.3）。在情景3中，所有響應變量具有固定的高共性（.7）。公共因子和唯一因子方差的確定方式與情景1類似。情景2和3中的參數生成分布 PΘ 是狄拉克測度（即，點質量集中在真實參數上）。我們使用 MATLAB 版本 23.2（MathWorks, 2023）生成數據；每種條件下運行了 512 次重復實驗。

抽樣與調參細節
我們采用了 Asparouhov 與 Muthén（2010）所推薦的先驗設定：對因子載荷使用非正常均勻先驗（improper uniform priors），對公共因子方差和唯一因子方差使用逆伽馬先驗 IG(1, 2)。MCMC 抽樣通過 JAGS（Plummer, 2017）完成。由于 JAGS 無法處理非正常先驗，我們改用一個彌散的正態先驗 N(0, 101?) 來代替因子載荷的先驗，該先驗與非正常均勻先驗幾乎無法區分。在每次重復實驗中，我們并行運行 5 條鏈。每條鏈的自適應迭代次數、預燒期（burn-in）迭代次數和保留迭代次數分別為 1000、10000 和 10000。利用這 5 條鏈的保留迭代樣本，以 10 的間隔進行抽稀（thinning），共獲得 5000 個 θ 的蒙特卡洛樣本。在每次重復中，記錄 θ 每個分量的潛在尺度縮減因子（PSRF）和有效樣本量（ESS）。若所有參數均滿足 PSRF ≤ 1.1 且 ESS ≥ 100，則認為該次重復實驗已收斂。

為評估后驗可能性并進行校準，我們計算了Wald檢驗統計量[11]和PDR統計量[13]。我們使用MATLAB內置的信賴域算法（通過fminunc）對關于θ的對數后驗進行最大化；我們提供了后驗的解析梯度和期望Hessian矩陣以加速優化過程。相同的期望Hessian矩陣也用于定義Wald檢驗統計量，以近似MAP估計量的協方差矩陣。

為調優SPRSA算法，我們進行了初步模擬實驗。具體而言，學習率序列由 a? = αk?? 確定，其中 α = 0.1，β = 0.65；有限差分率序列設定為 c? = γk??，其中 γ = 0.05，δ = 0.149。可以直觀驗證，這兩個速率序列滿足條件[22]。算法的迭代次數設定為 K = 50000。最后一次迭代的有限差分擾動值為 c????? = 0.05 × 50000??·1?? ≈ 0.01。校準后的α水平的平均估計值使用公式[23]計算得出。我們在MATLAB中實現了SPRSA算法；若本文被接受發表，我們將在補充材料中提供源代碼。

評估標準

我們評估了基于后驗推斷的貝葉斯有效性，無論是否經過校準。在每次重復實驗中，我們使用觀測數據 y 和數據生成參數 θ ~ P*Θ 來計算原始后驗等高線 ?y(θ) 和校準后的后驗等高線 ??y(θ)。隨后，我們在每種模擬條件下，針對原始后驗等高線值和校準后后驗等高線值，在所有重復實驗中分別構建了經驗分布函數（EDFs）。若EDF曲線位于對角線下方，則表明推斷是保守的，因而也是有效的；而若曲線位于對角線上方，則表明有效性要求未被滿足。為考慮蒙特卡洛誤差（MC error），與對角線的比較將參照其95%正態近似蒙特卡洛置信帶（即，

所有 θ 的聯合推斷結果總結于圖2中。當響應變量數量較小時（即，m = 5），基于橢圓和HPD可信區域的貝葉斯推斷往往無效；一個例外出現在情景3中的橢圓區域。當共性固定且較低時（情景2），橢圓區域表現得最為寬松；而當共性均勻生成時（情景1），HPD區域表現得最為寬松。相比之下，隨著響應變量數量增加至 m = 15，有效性通常得以恢復。但在情景3中，當共性固定且較高時，存在一個例外：基于HPD區域的推斷仍保持不可接受的寬松性。

同時，使用所提算法進行校準后的后驗可能性在所有模擬條件下均實現了有效性。然而，有效性保證會引入保守性，其程度取決于若干因素。特別是，校準后的貝葉斯推斷在響應變量較少（即，m = 5）和共性較低（即，情景2）時變得更加保守。校準后，HPD區域通常比橢圓區域更不保守。值得注意的是，當 m = 15 時，基于HPD區域的校準可能性等高線幾乎在所有三個參數生成情景下都達到了一致性。

對這一觀察的一個可能解釋是，與橢圓區域相比，HPD區域能更好地近似 πy 的上層水平集的形狀。

示例

為結束本節，我們說明校準結果在實踐中如何呈現。當分析單個數據集并選定檢驗統計量后，我們可以在一組預設的閾值 ξ?, ..., ξQ 上重復校準程序，從而得到一組校準后的 α 水平，α*(ξ?), ..., α*(ξQ)。為了突出校準的效果，我們建議選擇 {ξq} 作為在一組名義 α 水平下的檢驗統計量的 (1?α) 后驗分位數。然后可將校準后的 α 水平相對于后驗分位數的名義 α 水平作圖，創建一種類似于百分位-百分位圖的可視化圖形。

我們以在 m = 5 且參數按情景1生成條件下生成的最后一組數據為例。基于保留的MCMC抽樣，Wald統計量和PDR統計量的估計后驗密度顯示在圖3左側面板中。對于每個統計量，令 Q = 19，ξ?, ..., ξ?? 為后驗分布的 .95, .9, ..., .05 分位數。建議的圖形展示呈現在圖3右側面板中。在此數據集中，校準推斷對于兩種類型的可信區域而言，均比原始后驗推斷更為保守。更具體地說，對于HPD區域，α 水平調整幅度比橢圓區域更大。

貝葉斯統計在心理學家中廣受歡迎，因其能提供直觀的不確定性量化、適用于多種建模場景，并且在小樣本情況下有時表現優異（例如，Depaoli, 2021；Muthén & Asparouhov, 2012；van de Schoot 等, 2021）。然而，借助推斷模型（IM）中關于貝葉斯有效性的關鍵概念，我們表明：當用于推斷的先驗與真實參數生成機制（即真實的參數生成先驗）不匹配時，貝葉斯方法在重復抽樣下可能是不可靠的。由于在實際應用中，數據分析者通常無法獲知這一真實生成機制，因此我們提出一種更安全的替代方案：對基于后驗的推斷進行校準，以實現頻率學派的有效性（frequentist validity）——這是一種更強的要求，能夠保證在任意參數生成先驗下都滿足貝葉斯有效性。

為解決校準問題，我們開發了一種SPRSA算法，該算法將流形優化與無梯度隨機逼近（SA）相結合。隨后，我們報告了一項針對簡單單因子模型的蒙特卡洛實驗。結果表明，采用一種廣泛使用的先驗設定的標準貝葉斯推斷，其有效性會因可信區域類型、響應變量數量以及真實參數生成機制的不同而失效。相比之下，經過校準的貝葉斯推斷在所有模擬條件下均實現了有效性。此外，我們還證明了SPRSA算法可擴展至心理學應用中常見的現實問題規模。文中也提供了校準結果的建議可視化圖形展示方式。

本研究存在若干局限，有待未來研究加以解決。首先，我們的模擬僅限于一個簡單的單因子模型。鑒于貝葉斯方法的廣泛應用，有必要開展更全面的蒙特卡洛實驗，以評估常用先驗在多大程度上會導致無效推斷，并凸顯校準的普遍必要性。其次，Wald 統計量和 PDR 統計量的使用要求在 SPRSA 算法的每次迭代中求解兩次最大后驗（MAP）估計，對于復雜模型而言，這可能帶來較大的計算負擔。未來的研究可探索替代的有限差分（FD）梯度估計方法或檢驗統計量，以進一步減輕計算開銷。最后，我們的方法假設 p 值函數是可微的，因此無法直接應用于離散數據問題。一種有前景的解決方案是對檢驗統計量加入隨機擾動，從而強制其分布變為連續的。

https://www.researchgate.net/publication/397188720_Calibrating_Bayesian_Inference