全文連結:https://tecdat.cn/?p=33416
原文出處:拓端數據部落公眾號
介紹
在這裡,我們將幫助客戶將 PyMC3 用於兩個貝葉斯推理案例研究:拋硬幣和保險索賠發生。
方法:
回想一下,我們最初的貝葉斯推理方法是:
使用 PyMC3,我們現在可以簡化和壓縮這些步驟。
首先,我們設定先驗信念和先驗β-二項分布。
prior_beta = prior_beta.pdf(theta) / prior_beta.pdf(theta).sum() # 樣本積分 [pmf]()
plt.legend();
其次,我們定義並檢查我們的樣本觀察數據
print(f'Observed P(tails) = {tails/trials}')
第三,我們定義並運行我們的數學模型
請注意,PyMC3 提供了一種乾淨有效的語法來描述先驗分布和觀測數據,我們可以從中包括或單獨啟動模型抽樣。另請注意,PyMC3 允許我們定義先驗、引入樣本觀察數據並啟動後驗模擬。
# [NUTS](),採樣器(漢密爾頓式)
step = pm.NUTS()
結果
或者通過更多的採樣和更多的鏈。然後,跟蹤摘要返回有用的模型性能摘要統計信息:
- mc_error通過將跡線分解為批次,計算每個批次的平均值,然後計算這些平均值的標準偏差來估計模擬誤差。
- hpd_* 給出最高的後密度區間。2.5 和 97.5 標籤有點誤導。有很多 95% 的可信區間,具體取決於左右尾巴的相對權重。95% HPD 區間是這 95% 區間中最窄的。
- Rhat有時被稱為潛在的規模縮減因子,它為我們提供了一個因子,如果我們的MCMC鏈更長,則可以減少方差。它是根據鏈與每個鏈內的方差來計算的。接近 1 的值很好。
summary
我們使用跡線手動繪製和比較先驗分布和後驗分布。確認這些與手動獲得的相似,後驗分布均值為 P(Tails|觀測數據)= 0.35。
但是,PyMC3還提供了創建跡線圖,後驗分布圖。
pm.traceplot(trace)
pm.plot_posterior(trace,ref_val=0.5);
我們有它。PyMC3 和其他類似軟體包提供了一組簡單的函數來組裝和運行機率模擬,例如貝葉斯推理。
個案研究:
使用貝葉斯推理評估保險索賠發生率
保險索賠通常被建模為由於泊松分布式過程而發生。
泊松分布由下式給出:
其中 lambda λ 是事件的「速率」,由事件總數 (k) 除以數據中的單位數 (n) 給出 (λ = k/n)。在泊松分布中,泊松分布的期望值 E(Y)、均值 E(X) 和方差 Var(Y) 相同;
例如,E(Y) = E(X) = Var(X) = λ。
請注意,如果方差大於均值,則稱數據過於分散。這在具有大量零的保險索賠數據中很常見,並且最好由負二項式和零膨脹模型(如 ZIP 和 ZINB)處理。
一、建立先驗分布
在這裡,我們生成一些觀測數據,這些數據遵循泊松分布,速率為 lambda,λ = 2。
n = 1000
lam_ = 2
axs.set_title('Histogram: Simulated Poisson $y$')
axs.set_xlabel('Poisson lambda=λ')
axs.set_ylabel('P(λ)')
axs.legend();
我們可以使用β泊松,或任何類似於觀察到的λ數據形狀的分布,但是伽馬泊松最適合:
- 泊松可以取任何正數到無窮大(0,∞),而β或均勻是[0-100]。
- 伽馬和泊松屬於同一分布家族。
- 伽馬的峰值接近於零。
- 伽馬尾巴走向無窮大。
伽馬泊松先驗為:
其中 a 是伽馬形狀,b 是伽馬速率參數。伽馬密度函數為:
其中 a>0 是形狀參數,b>0 是速率參數,以及
和
注意在 scipy 中,伽馬分布使用形狀 a 和尺度參數化,其中速率 b 等於尺度的倒數(速率 = 1/尺度)。
prior = lambda x: stats.gamma.pdf(x, a=a, scale=rate,loc=0)
priors = prior(x)
# 畫圖
axs.plot(x, priors, 'r-',label='Gamma')
二、似然函數與後驗
伽馬函數通常被稱為廣義階乘,因為:
sp.gamma(n+1) == math.factorial(n)
True
則似然函數為:
然後作為
後向分布再次為伽馬
def posterior(lam,y):
shape = a + y.sum()
如圖所示,後驗平均值(藍色)以我們在開始時設置的真實 lambda 速率為中心。後驗平均值為:
即後驗平均值是先驗平均值和觀測樣本平均值的加權平均值
posterior mean: {(a+y.sum()) / (b+y.size)}
sample mean:{y.mean()}""")
現在讓我們在 PyMC3 中重現上述步驟。
print(a,b,lam_,y.shape)
with model:
# 定義參數 [lambda]() 的先驗值。
prior_lam = pm.Gamma('prior-gamma-lambda', alpha=a, beta=b)
跡線圖顯示每個模擬的結果。
低於平均值、分位數、可信區間 (HPD) 94% 和任意參考值(橙色垂直)。
import warnings
with warnings.catch_warnings():
warnings.simplefilter("ignore")
您可能已經注意到,在這個例子中,我們已經根據觀察到的數據定義了我們的先驗分布,並對該數據應用貝葉斯推理來推導出後驗分布,確認 lambda 為 2。
結論:
在這篇文章中,PyMC3 被應用於對兩個示例進行貝葉斯推理:使用 β-二項分布的拋硬幣偏差,以及使用 gamma-泊松分布的保險索賠發生。
最受歡迎的見解
1.matlab使用貝葉斯優化的深度學習
2.matlab貝葉斯隱馬爾可夫hmm模型實現
3.R語言Gibbs抽樣的貝葉斯簡單線性回歸仿真
4.R語言中的block Gibbs吉布斯採樣貝葉斯多元線性回歸
5.R語言中的Stan機率編程MCMC採樣的貝葉斯模型
6.R語言貝葉斯Poisson泊松-正態分布模型分析職業足球比賽進球數
7.R語言使用貝葉斯 層次模型進行空間數據分析
8.R語言隨機搜索變量選擇SSVS估計貝葉斯向量自回歸(BVAR)模型
9.matlab貝葉斯隱馬爾可夫hmm模型實現