全文鏈接：https://tecdat.cn/?p=33416

原文出處：拓端數(shù)據(jù)部落公眾號(hào)

介紹

在這里，我們將幫助客戶將 PyMC3 用于兩個(gè)貝葉斯推理案例研究：拋硬幣和保險(xiǎn)索賠發(fā)生。

方法：

回想一下，我們最初的貝葉斯推理方法是：

1. 設(shè)置先前的假設(shè)，并根據(jù)啟發(fā)式、歷史或樣本數(shù)據(jù)建立我們數(shù)據(jù)的“已知已知”。

2. 形式化問(wèn)題空間和先前假設(shè)的數(shù)學(xué)模型。

3. 正式化先前的分布。

4. 應(yīng)用貝葉定理從觀察到的樣本數(shù)據(jù)中推導(dǎo)出后驗(yàn)參數(shù)值。

5. 重復(fù)步驟 1-4，以獲取更多數(shù)據(jù)樣本。

使用 PyMC3，我們現(xiàn)在可以簡(jiǎn)化和壓縮這些步驟。

首先，我們?cè)O(shè)定先驗(yàn)信念和先驗(yàn)β-二項(xiàng)分布。

prior_beta = prior_beta.pdf(theta) / prior_beta.pdf(theta).sum() # 樣本積分?[pmf]()plt.legend();

其次，我們定義并檢查我們的樣本觀察數(shù)據(jù)

print(f'Observed P(tails) = {tails/trials}')

第三，我們定義并運(yùn)行我們的數(shù)學(xué)模型

請(qǐng)注意，PyMC3 提供了一種干凈有效的語(yǔ)法來(lái)描述先驗(yàn)分布和觀測(cè)數(shù)據(jù)，我們可以從中包括或單獨(dú)啟動(dòng)模型抽樣。另請(qǐng)注意，PyMC3 允許我們定義先驗(yàn)、引入樣本觀察數(shù)據(jù)并啟動(dòng)后驗(yàn)?zāi)M。

? ? ? ?# [NUTS]()，采樣器（漢密爾頓式） ? ?step = pm.NUTS() ? ?

結(jié)果

或者通過(guò)更多的采樣和更多的鏈。然后，跟蹤摘要返回有用的模型性能摘要統(tǒng)計(jì)信息：

• mc_error通過(guò)將跡線分解為批次，計(jì)算每個(gè)批次的平均值，然后計(jì)算這些平均值的標(biāo)準(zhǔn)偏差來(lái)估計(jì)模擬誤差。

• hpd_* 給出最高的后密度區(qū)間。2.5 和 97.5 標(biāo)簽有點(diǎn)誤導(dǎo)。有很多 95% 的可信區(qū)間，具體取決于左右尾巴的相對(duì)權(quán)重。95% HPD 區(qū)間是這 95% 區(qū)間中最窄的。

• Rhat有時(shí)被稱為潛在的規(guī)?？s減因子，它為我們提供了一個(gè)因子，如果我們的MCMC鏈更長(zhǎng)，則可以減少方差。它是根據(jù)鏈與每個(gè)鏈內(nèi)的方差來(lái)計(jì)算的。接近 1 的值很好。

summary

我們使用跡線手動(dòng)繪制和比較先驗(yàn)分布和后驗(yàn)分布。確認(rèn)這些與手動(dòng)獲得的相似，后驗(yàn)分布均值為 P（Tails|觀測(cè)數(shù)據(jù)）= 0.35。

但是，PyMC3還提供了創(chuàng)建跡線圖，后驗(yàn)分布圖。

? ?pm.traceplot(trace) ? ?pm.plot_posterior(trace,ref_val=0.5);

我們有它。PyMC3 和其他類似軟件包提供了一組簡(jiǎn)單的函數(shù)來(lái)組裝和運(yùn)行概率模擬，例如貝葉斯推理。

個(gè)案研究：

使用貝葉斯推理評(píng)估保險(xiǎn)索賠發(fā)生率

保險(xiǎn)索賠通常被建模為由于泊松分布式過(guò)程而發(fā)生。

泊松分布由下式給出：

其中 lambda λ 是事件的“速率”，由事件總數(shù) （k） 除以數(shù)據(jù)中的單位數(shù) （n） 給出 （λ = k/n）。在泊松分布中，泊松分布的期望值 E（Y）、均值 E（X） 和方差 Var（Y） 相同;

例如，E（Y） = E（X） = Var（X） = λ。

請(qǐng)注意，如果方差大于均值，則稱數(shù)據(jù)過(guò)于分散。這在具有大量零的保險(xiǎn)索賠數(shù)據(jù)中很常見(jiàn)，并且最好由負(fù)二項(xiàng)式和零膨脹模型（如 ZIP 和 ZINB）處理。

一、建立先驗(yàn)分布

在這里，我們生成一些觀測(cè)數(shù)據(jù)，這些數(shù)據(jù)遵循泊松分布，速率為 lambda，λ = 2。

n = 1000lam_ = 2axs.set_title('Histogram: Simulated Poisson $y$')axs.set_xlabel('Poisson lambda=λ')axs.set_ylabel('P(λ)')axs.legend();

我們可以使用β泊松，或任何類似于觀察到的λ數(shù)據(jù)形狀的分布，但是伽馬泊松最適合：

• 泊松可以取任何正數(shù)到無(wú)窮大（0，∞），而β或均勻是[0-100]。

• 伽馬和泊松屬于同一分布家族。

• 伽馬的峰值接近于零。

• 伽馬尾巴走向無(wú)窮大。

伽馬泊松先驗(yàn)為：

其中 a 是伽馬形狀，b 是伽馬速率參數(shù)。伽馬密度函數(shù)為：

其中 a>0 是形狀參數(shù)，b>0 是速率參數(shù)，以及

和

注意在 scipy 中，伽馬分布使用形狀?a?和尺度參數(shù)化，其中速率?b?等于尺度的倒數(shù)（速率 = 1/尺度）。

prior = lambda x: stats.gamma.pdf(x, a=a, scale=rate,loc=0)priors = prior(x)# 畫圖axs.plot(x, priors, 'r-',label='Gamma')

二、似然函數(shù)與后驗(yàn)

伽馬函數(shù)通常被稱為廣義階乘，因?yàn)椋?/p>

sp.gamma(n+1) == math.factorial(n)

True

則似然函數(shù)為：

然后作為

后向分布再次為伽馬?

def posterior(lam,y): ? ? ? ?shape = a + y.sum()

如圖所示，后驗(yàn)平均值（藍(lán)色）以我們?cè)陂_(kāi)始時(shí)設(shè)置的真實(shí) lambda 速率為中心。后驗(yàn)平均值為：

即后驗(yàn)平均值是先驗(yàn)平均值和觀測(cè)樣本平均值的加權(quán)平均值

posterior mean: {(a+y.sum()) / (b+y.size)}sample mean:{y.mean()}""")

現(xiàn)在讓我們?cè)?PyMC3 中重現(xiàn)上述步驟。

print(a,b,lam_,y.shape)

with model: ? ? ? ?# 定義參數(shù)?[lambda]()?的先驗(yàn)值。 ? ?prior_lam = pm.Gamma('prior-gamma-lambda', alpha=a, beta=b)

跡線圖顯示每個(gè)模擬的結(jié)果。

低于平均值、分位數(shù)、可信區(qū)間 （HPD） 94% 和任意參考值（橙色垂直）。

import warningswith warnings.catch_warnings(): ? ?warnings.simplefilter("ignore")

您可能已經(jīng)注意到，在這個(gè)例子中，我們已經(jīng)根據(jù)觀察到的數(shù)據(jù)定義了我們的先驗(yàn)分布，并對(duì)該數(shù)據(jù)應(yīng)用貝葉斯推理來(lái)推導(dǎo)出后驗(yàn)分布，確認(rèn) lambda 為 2。

結(jié)論：

在這篇文章中，PyMC3 被應(yīng)用于對(duì)兩個(gè)示例進(jìn)行貝葉斯推理：使用 β-二項(xiàng)分布的拋硬幣偏差，以及使用 gamma-泊松分布的保險(xiǎn)索賠發(fā)生。

?最受歡迎的見(jiàn)解

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)