Python用PyMC3貝葉斯模型平均BMA：采樣、信息準(zhǔn)則比較和預(yù)測(cè)可視化靈長(zhǎng)類動(dòng)物的乳汁成
全文鏈接：https://tecdat.cn/?p=33449
原文出處：拓端數(shù)據(jù)部落公眾號(hào)
當(dāng)面對(duì)多個(gè)模型時(shí)，我們有多種選擇。模型選擇因其簡(jiǎn)單性而具有吸引力，但我們正在丟棄有關(guān)模型中不確定性的信息。
print(f"Runing
模型平均
一種替代方法是執(zhí)行模型選擇，但討論所有不同的模型以及給定信息準(zhǔn)則的計(jì)算值。重要的是要將所有這些數(shù)字和測(cè)試放在我們問(wèn)題的背景下，以便我們和我們的客戶能夠更好地了解我們方法可能存在的局限性和缺點(diǎn)。如果你在學(xué)術(shù)界，你可以使用這種方法向論文、演示文稿、論文等的討論部分添加元素。
另一種方法是執(zhí)行模型平均?，F(xiàn)在的想法是使用模型的加權(quán)平均值生成元模型（和元預(yù)測(cè)）。有幾種方法可以做到這一點(diǎn)，PyMC3 包括其中的 3 種，我們將簡(jiǎn)要討論，您將在?Yuling Yao 等人的工作中找到更徹底的解釋。
偽貝葉斯模型平均
貝葉斯模型可以通過(guò)其邊緣概率進(jìn)行加權(quán)，這被稱為貝葉斯模型平均。我們可以使用以下公式來(lái)做到這一點(diǎn)：
這種方法稱為偽貝葉斯模型平均或類似赤池的加權(quán)，是一種啟發(fā)式方法，用于根據(jù)信息標(biāo)準(zhǔn)值計(jì)算每個(gè)模型（給定一組固定的模型）的相對(duì)概率?？纯捶帜钢皇且粋€(gè)歸一化項(xiàng)，以確保權(quán)重總和為 1。
使用貝葉斯自舉進(jìn)行偽貝葉斯模型平均
上述計(jì)算權(quán)重的公式是一種非常好且簡(jiǎn)單的方法，但有一個(gè)主要警告，它沒(méi)有考慮 IC 計(jì)算中的不確定性。
堆疊
在PyMC3中實(shí)現(xiàn)的第三種方法被稱為預(yù)測(cè)分布的堆疊，并且最近被提出。我們希望在一個(gè)元模型中組合多個(gè)模型，以最小化元模型和真實(shí)生成模型之間的分歧，當(dāng)使用對(duì)數(shù)評(píng)分規(guī)則時(shí)，這相當(dāng)于：
加權(quán)后驗(yàn)預(yù)測(cè)樣本
一旦我們計(jì)算了權(quán)重，使用上述 3 種方法中的任何一種，我們就可以使用它們來(lái)獲得加權(quán)后驗(yàn)預(yù)測(cè)樣本。PyMC3 提供了以簡(jiǎn)單方式執(zhí)行這些步驟的函數(shù)，因此讓我們通過(guò)示例查看它們的實(shí)際效果。
簡(jiǎn)而言之，我們的問(wèn)題如下：我們想探索幾種靈長(zhǎng)類動(dòng)物的乳汁成分，假設(shè)來(lái)自大腦較大的靈長(zhǎng)類動(dòng)物的雌性產(chǎn)生更有營(yíng)養(yǎng)的牛奶（這樣做是為了*支持這種大大腦的發(fā)育）。對(duì)于進(jìn)化生物學(xué)家來(lái)說(shuō)，這是一個(gè)重要的問(wèn)題，為了給出和回答，我們將使用3個(gè)變量，兩個(gè)預(yù)測(cè)變量：新皮層的比例與總質(zhì)量的比較 大腦和母親體重的對(duì)數(shù)。對(duì)于預(yù)測(cè)變量，每克牛奶的千卡。使用這些變量，我們將構(gòu)建 3 個(gè)不同的線性模型：
僅使用新皮層變量的模型
僅使用質(zhì)量變量對(duì)數(shù)的模型
使用兩個(gè)變量的模型
d.iloc[:, 1:] = d.iloc[:, 1:] - d.iloc[:, 1:].mean()d.head()
現(xiàn)在我們有了數(shù)據(jù)，我們將僅使用?neocortex。
with pm.Model() as model_0: ? ? ?trace_0 = pm.sample(2000, return_inferencedata=True)
第二個(gè)模型與第一個(gè)模型完全相同，只是我們現(xiàn)在使用質(zhì)量的對(duì)數(shù)
with pm.Model() as model_1: ? ?trace_1 = pm.sample(2000, return_inferencedata=True)
最后是第三個(gè)模型使用?neocortex和 變量log_mass
with pm.Model() as model_2: ? ? ?trace_2 = pm.sample(2000, return_inferencedata=True)
現(xiàn)在我們已經(jīng)對(duì) 3 個(gè)模型的后驗(yàn)進(jìn)行了采樣，我們將對(duì)它們進(jìn)行視覺比較。一種選擇是使用forestplot支持繪制多個(gè)跡線的函數(shù)。
az.plot_fo
?另一種選擇是在同一圖中繪制多條跡線是使用densityplot?。
az.plot_d
?現(xiàn)在我們已經(jīng)對(duì) 3 個(gè)模型的后驗(yàn)進(jìn)行了采樣，我們將使用 WAIC（廣泛適用的信息標(biāo)準(zhǔn)）來(lái)比較 3 個(gè)模型。我們可以使用 PyMC3 附帶的compare功能來(lái)做到這一點(diǎn)。
comp = az.compare(model_dict)comp
我們可以看到最好的模型是，具有兩個(gè)預(yù)測(cè)變量的模型。請(qǐng)注意，數(shù)據(jù)幀按從最低到最高 WAIC 的順序（即從好到最差的模型）。
現(xiàn)在，我們將使用copmuted來(lái)生成預(yù)測(cè)，而不是基于單個(gè)模型，而是基于加權(quán)模型集。
ppc_w = pm.sample_posterior_predictive_w( ?
請(qǐng)注意，我們正在傳遞按其索引排序的權(quán)重。
我們還將計(jì)算最低 WAIC 模型的 PPC
ppc_2 = pm.sample_posterior_predi
比較這兩種預(yù)測(cè)的一種簡(jiǎn)單方法是繪制它們的平均值和 hpd 區(qū)間
plt.yticks([])plt.ylim(-1, 2)plt.legend();
正如我們所看到的，兩個(gè)預(yù)測(cè)的平均值幾乎相同，但加權(quán)模型中的不確定性更大。我們已經(jīng)有效地將我們應(yīng)該選擇哪個(gè)模型的不確定性傳遞到后驗(yàn)預(yù)測(cè)樣本中。
結(jié)語(yǔ)：
還有其他方法可以平均模型，例如，顯式構(gòu)建一個(gè)包含我們擁有的所有模型的元模型。然后，我們?cè)谀Ｐ椭g跳轉(zhuǎn)時(shí)執(zhí)行參數(shù)推理。這種方法的一個(gè)問(wèn)題是，在模型之間跳躍可能會(huì)妨礙后驗(yàn)的正確采樣。
版本信息
%load_ext watermark%watermark -n -u -v -iv -w
最受歡迎的見解
1.matlab使用貝葉斯優(yōu)化的深度學(xué)習(xí)
2.matlab貝葉斯隱馬爾可夫hmm模型實(shí)現(xiàn)
3.R語(yǔ)言Gibbs抽樣的貝葉斯簡(jiǎn)單線性回歸仿真
4.R語(yǔ)言中的block Gibbs吉布斯采樣貝葉斯多元線性回歸
5.R語(yǔ)言中的Stan概率編程MCMC采樣的貝葉斯模型
6.R語(yǔ)言貝葉斯Poisson泊松-正態(tài)分布模型分析職業(yè)足球比賽進(jìn)球數(shù)
7.R語(yǔ)言使用貝葉斯 層次模型進(jìn)行空間數(shù)據(jù)分析
8.R語(yǔ)言隨機(jī)搜索變量選擇SSVS估計(jì)貝葉斯向量自回歸（BVAR）模型
9.matlab貝葉斯隱馬爾可夫hmm模型實(shí)現(xiàn)
