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

2023-08-15 tecdat拓端

全文連結：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");

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

總結