Python用PyMC貝葉斯GLM廣義線性模型、NUTS採樣器擬合、後驗分布可視化

2023-08-15     tecdat拓端

原標題:Python用PyMC貝葉斯GLM廣義線性模型、NUTS採樣器擬合、後驗分布可視化

全文連結:https://tecdat.cn/?p=33436

原文出處:拓端數據部落公眾號

儘管貝葉斯方法相對於頻率主義方法的理論優勢已經在其他地方進行了詳細討論,但其更廣泛採用的主要障礙是「可用性」。而使用貝葉斯方法,客戶可以按照自己認為合適的方式定義模型。

線性回歸

在此示例中,我們將幫助客戶從最簡單的 GLM – 線性回歸開始。 一般來說,頻率論者對線性回歸的看法如下:

然後,我們可以使用普通最小二乘法(OLS)或最大似然法來找到最佳擬合。

機率重構

貝葉斯主義者對世界採取機率觀,並用機率分布來表達這個模型。我們上面的線性回歸可以重新表述為:

換句話說,我們將Y其視為一個隨機變量(或隨機向量),其中每個元素(數據點)都根據正態分布分布。此正態分布的均值由具有方差sigma的線性預測變量提供。

PyMC 中的貝葉斯 GLM

要開始在 PyMC 中構建 GLM,讓我們首先導入所需的模塊。

print(f"Running on PyMC v{pm.__version__}")

az.style.use("arviz-darkgrid")

數據

本質上,我們正在創建一條由截距和斜率定義的回歸線,並通過從均值設置為回歸線的正態採樣來添加數據點。

y = true_regression_line + rng.normal(scale=0.5, size=size)

data = pd.DataFrame(dict(x=x, y=y))

plt.legend(loc=0);

估計模型

讓我們將貝葉斯線性回歸模型擬合到此數據。

# 定義似然函數

likelihood = Normal("y", mu=intercept + slope * x, sigma=sigma, observed=y)

# 使用NUTS採樣推斷

idata = sample(3000)

對於了解機率編程的人來說,這應該是相當可讀的。

import bambi as bmb

idata = model.fit(draws=3000)

要短得多,但這段代碼與之前的規範完全相同(如果我們願意,您也可以更改先驗和其他所有內容)。

分析模型

貝葉斯推理不僅給了我們一條最佳擬合線(就像最大似然那樣),而是給出了合理參數的整個後驗分布。讓我們繪製參數的後驗分布和我們繪製的單個樣本。

az.plot_trace(idata, figsize=(10, 7));

左側顯示了我們的邊緣後驗 – 對於 x 軸上的每個參數值,我們在 y 軸上得到一個機率,告訴我們該參數值的可能性。

首先,各個參數(左側)的採樣鏈看起來均勻且平穩(沒有大的漂移或其他奇怪的模式)。

其次,每個變量的最大後驗估計值(左側分布中的峰值)非常接近用於生成數據的真實參數(x是回歸係數,sigma是我們正態的標準差)。

因此,在 GLM 中,我們不僅有一條最佳擬合回歸線,而且有許多。後驗預測圖從後驗圖(截距和斜率)中獲取多個樣本,並為每個樣本繪製一條回歸線。我們可以直接使用後驗樣本手動生成這些回歸線。

idata.posterior["y_model"] = idata.posterior["Intercept"] + idata.posterior["x"] * xr.DataArray(x)

_, ax = plt.subplots(figsize=(7, 7))

az.plot_lm(idata=idata, y="y", num_samples=100, axes=ax, y_model="y_model")

ax.set_title("Posterior predictive regression lines")

ax.set_xlabel("x");

我們估計的回歸線與真正的回歸線非常相似。但是由於我們只有有限的數據,我們的估計存在不確定性,這裡用線的可變性來表示。

總結

  • 可用性目前是更廣泛採用貝葉斯統計的巨大障礙。
  • Bambi允許使用從 R 借用的便捷語法進行 GLM 規範。然後可以使用pymc 進行推理。
  • 後驗預測圖使我們能夠評估擬合度和其中的不確定性。

延伸閱讀

有關其他背景信息,以下是一些關於貝葉斯統計的好資源:

  • 約翰·克魯施克(John Kruschke)的優秀著作《做貝葉斯數據分析》。

版本信息:

%load_ext watermark

%watermark -n -u -v -iv -w -p pytensor

Python implementation: CPython

Python version : 3.11.4

IPython version : 8.14.0

pytensor: 2.14.2

pymc : 5.7.2+0.gd59a960f.dirty

bambi : 0.12.0

arviz : 0.16.1

xarray : 2023.7.0

matplotlib: 3.7.2

numpy : 1.25.2

sys : 3.11.4 | packaged by conda-forge | (main, Jun 10 2023, 18:08:17) [GCC 12.2.0]

pandas : 2.0.3

Watermark: 2.4.3

最受歡迎的見解

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模型實現

文章來源: https://twgreatdaily.com/zh/58948dd27e9f82939f29ca95f38e467d.html