在PyStan中應用貝葉斯回歸

鼓勵您在參與本文之前檢查一下此概念性背景。

設定

Stan [1]是用于貝葉斯模型擬合的計算引擎。 它依賴于哈密頓量的蒙特卡羅(HMC)[2]的變體來從各種貝葉斯模型的后驗分布中采樣。

以下是設置Stan的詳細安裝步驟：

對于MacOS：

· 安裝miniconda / anaconda

· 安裝xcode

· 更新您的C編譯器：conda install clang_osx-64 clangxx_osx-64 -c anaconda

· 創建一個環境stan或pystan

· 輸入conda install numpy來安裝numpy或將numpy替換為您需要安裝的軟件包

· 安裝pystan：conda install -c conda-forge pystan

· 或者：pip install pystan

· 同時安裝：arviz，matplotlib，pandas，scipy，seaborn，statsmodels，pickle，scikit-learn，nb_conda和nb_conda_kernels

設置完成后，我們可以打開(Jupyter)筆記本并開始我們的工作。 首先，讓我們使用以下內容導入庫：

import pystanimport pickleimport numpy as npimport arviz as azimport pandas as pdimport seaborn as snsimport statsmodels.api as statmodimport matplotlib.pyplot as pltfrom IPython.display import Imagefrom IPython.core.display import HTML

投幣推理

回想一下在《概念導論》中我談到過如何在地面上找到一枚硬幣并將其拋擲K次以檢查它是否公平。

我們選擇了以下模型：

下面的概率質量函數對應于Y：

首先，通過使用NumPy的隨機功能進行仿真，我們可以了解參數為a = b = 5的先驗行為：

sns.distplot(np.random.beta(5,5, size=10000),kde=False)

如概念背景中所述，此先驗似乎是合理的，因為其對稱性約為0.5，但仍可能在任一方向上產生偏差。 然后，我們可以使用以下語法在pystan中定義模型：

· 數據對應于我們模型的數據。 在這種情況下，整數N對應于拋硬幣的次數，y對應于長度為N的整數的向量，該向量將包含我的實驗觀察結果。

· 參數對應于我們模型的參數，在這種情況下為theta或獲得"頭部"的概率。

· 模型對應于我們的先驗(β)和可能性(bernoulli)的定義。

# bernoulli modelmodel_code = """ data { int N; int y[N]; } parameters { real theta; } model { theta ~ beta(5, 5); for (n in 1:N) y[n] ~ bernoulli(theta); } """data = dict(N=4, y=[0, 0, 0, 0])model = pystan.StanModel(model_code=model_code)fit = model.sampling(data=data,iter=4000, chains=4, warmup=1000)la = fit.extract(permuted=True) # return a dictionary of arraysprint(fit.stansummary())

請注意，model.sampling的默認參數是iter = 1000，chains = 4和warmup = 500。 我們會根據時間和可用的計算資源進行調整。

· iter≥0對應于我們的每條MCMC鏈的運行次數(對于大多數應用程序，其運行次數不得少于1,000)

· warmup或"burn-in"≥0對應于我們采樣開始時的初始運行次數。 考慮到這些鏈條在運行之初就非常不穩定和幼稚，實際上，我們通常將這個量定義為丟棄第一個B =warmup樣本數，否則我們會在估計中引入不必要的噪聲。

· chains≥0對應于我們樣本中的MCMC鏈數。

以下是上述模型的輸出：

· 我們θ的后均值約為0.36 <0.5。 盡管如此，95％的后可信區間很寬：(0.14，0.61)。 因此，我們可以認為結果在統計上不是結論性的，但它指向偏見而沒有跳到0。

應用回歸：每加侖汽車和英里數(MPG)

Image by Susan Sewert from Pixabay

讓我們構建一個貝葉斯線性回歸模型來解釋和預測汽車數據集中的MPG。 盡管我的方法是傳統的基于可能性的模型，但我們可以(并且應該！)使用原始ML訓練檢驗分裂來評估我們的預測質量。

該數據集來自UCI ML存儲庫，包含以下信息：

我們為什么不堅持我們的基礎知識并使用標準的線性回歸？ [3]回顧其以下功能形式：

現在，對于貝葉斯公式：

盡管前兩行看起來完全一樣，但是現在我們需要為β和σ(以及α，為了表示簡單起見，我選擇將其吸收到β向量中)建立先驗分布。

您可以在此處詢問有關Σ的信息。 這是一個N(訓練觀察數)x P(系數/特征數)矩陣。 在標準線性回歸情況下，要求此Σ是對角矩陣，且沿對角線有σ(觀測值之間的獨立性)。

在執行EDA時，您應該始終考慮以下幾點：

· 我的可能性是多少？

· 我的模特應該是什么？ 互動與沒有互動？ 多項式？

· 我的參數(和超參數)是什么？應該選擇哪種先驗條件？

· 我應該考慮任何聚類，時間或空間依賴性嗎？

EDA

與任何類型的數據分析一樣，至關重要的是首先了解我們感興趣的變量/功能以及它們之間的相互關系。

首先加載數據集：

cars_data = pd.read_csv("~/cars.csv").set_index("name")print(cars_data.shape)cars_data.head()

讓我們檢查一下目標變量和預測變量之間的關系(如果有)：

There appears to be some interesting separation between American and Japanese cars w.r.t. mpg.

現在，讓我們檢查一下數字變量和mpg之間的關系：

· 我更喜歡形象化這些關系而不是計算相關性，因為它給了我關于它們之間關系的視覺和數學意義，超越了簡單的標量。

All but year and acceleration are negatively related to mpg, i.e., for an increase of 1 unit in displacement, there appears to be a decrease in mpg.

訓練/擬合

讓我們準備數據集以進行擬合和測試：

from numpy import randomfrom sklearn import preprocessing, metrics, linear_modelfrom sklearn.model_selection import train_test_splitfrom sklearn.metrics import mean_squared_errorrandom.seed(12345)cars_data = cars_data.set_index('name')y = cars_data['mpg']X = cars_data.loc[:, cars_data.columns != 'mpg']X = X.loc[:, X.columns != 'name']X = pd.get_dummies(X, prefix_sep='_', drop_first=False) X = X.drop(columns=["origin_European"]) # This is our reference categoryX_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.15, random_state=0)X_train.head()

現在進入我們的模型規范，它與上面的拋硬幣問題沒有實質性的區別，除了我使用矩陣符號來盡可能簡化模型表達式的事實之外：

# Succinct matrix notationcars_code = """data { int N; // number of training samples int K; // number of predictors - 1 (intercept) matrix[N, K] x; // matrix of predictors vector[N] y_obs; // observed/training mpg int N_new; matrix[N_new, K] x_new;}parameters { real alpha; vector[K] beta; real sigma; vector[N_new] y_new;}transformed parameters { vector[N] theta; theta = alpha + x * beta;}model { sigma ~ exponential(1); alpha ~ normal(0, 6); beta ~ multi_normal(rep_vector(0, K), diag_matrix(rep_vector(1, K))); y_obs ~ normal(theta, sigma); y_new ~ normal(alpha + x_new * beta, sigma); // prediction model}"""

· 數據與我們模型的數據相對應，如上面的拋硬幣示例所示。

· 參數對應于我們模型的參數。

· 此處經過轉換的參數使我可以將theta定義為模型在訓練集上的擬合值

· 模型對應于我們對sigma，alpha和beta的先驗定義，以及我們對P(Y | X，α，β，σ)的似然性(正常)的定義。

在對初始數據進行檢查之后，選擇了以上先驗條件：

· 為什么在σ上有指數先驗？ 好吧，σ≥0(根據定義)。 為什么不穿制服或伽瑪呢？ PC框架[4]-我的目標是最簡約的模型。

· 那么α和β呢？ 在正常情況下，最簡單的方法是為這些參數選擇常規先驗。

· 為什么要對β使用multi_normal()？ 數學統計中有一個眾所周知的結果，該結果表明，長度為

stan的有趣之處在于，我們可以要求在測試集上進行預測而無需重新擬合-而是可以在第一個調用中從stan請求它們。

· 我們通過在上面的單元格中定義t_new來實現這一點。

· theta是我們對訓練集的合適預測。

我們指定要提取的數據集，并繼續從模型中取樣，同時將其保存以備將來參考：

cars_dat = {'N': X_train.shape[0], 'N_new': X_test.shape[0], 'K': X_train.shape[1], 'y_obs': y_train.values.tolist(), 'x': np.array(X_train), 'x_new': np.array(X_test)}sm = pystan.StanModel(model_code=cars_code)fit = sm.sampling(data=cars_dat, iter=6000, chains=8)# Save fitted model!with open('bayes-cars.pkl', 'wb') as f: pickle.dump(sm, f, protocol=pickle.HIGHEST_PROTOCOL)# Extract and print the output of our modella = fit.extract(permuted=True)print(fit.stansummary())

這是我打印的輸出(出于本教程的目的，我擺脫了That和n_eff)：

Inference for Stan model: anon_model_3112a6cce1c41eead6e39aa4b53ccc8b.8 chains, each with iter=6000; warmup=3000; thin=1; post-warmup draws per chain=3000, total post-warmup draws=24000.mean se_mean sd 2.5% 25% 50% 75% 97.5% alpha -8.35 0.02 3.81 -15.79 -10.93 -8.37 -5.74 -0.97 beta[1] -0.48 1.9e-3 0.33 -1.12 -0.7 -0.48 -0.26 0.17 beta[2] 0.02 4.6e-5 7.9e-3 6.1e-3 0.02 0.02 0.03 0.04 beta[3] -0.02 8.7e-5 0.01 -0.04 -0.03 -0.02-6.9e-3 0.01 beta[4] -7.0e-3 4.0e-6 7.0e-4-8.3e-3-7.4e-3-7.0e-3-6.5e-3-5.6e-3 beta[5] 0.1 5.9e-4 0.1 -0.1 0.03 0.1 0.17 0.3 beta[6] 0.69 2.7e-4 0.05 0.6 0.65 0.69 0.72 0.78 beta[7] -1.75 2.9e-3 0.51 -2.73 -2.09 -1.75 -1.41 -0.73 beta[8] 0.36 2.7e-3 0.51 -0.64 0.02 0.36 0.71 1.38 sigma 3.33 7.6e-4 0.13 3.09 3.24 3.32 3.41 3.59 y_new[1] 26.13 0.02 3.37 19.59 23.85 26.11 28.39 32.75 y_new[2] 25.38 0.02 3.36 18.83 23.12 25.37 27.65 31.95......theta[1] 24.75 2.7e-3 0.47 23.83 24.44 24.75 25.07 25.68 theta[2] 19.59 2.6e-3 0.43 18.76 19.3 19.59 19.88 20.43 ......

fit.stansummary()是一個像表一樣排列的字符串，它為我提供了擬合期間估計的每個參數的后驗均值，標準差和幾個百分點。 雖然alpha對應于截距，但我們有8個beta回歸變量，復雜度或觀察誤差sigma，訓練集theta的擬合值以及測試集y_new的預測值。

診斷

在幕后運行MCMC，至關重要的是我們要以圖表的形式檢查基本診斷。 對于這些情節，我依靠出色的arviz庫進行貝葉斯可視化(與PyMC3和pystan同樣有效)。

· 鏈混合-軌跡圖。 這些系列圖應顯示在實線中"合理大小"的間隔內振蕩的"厚發"鏈，以指示良好的混合。 "稀疏的"鏈意味著MCMC無法有效地進行探索，并且可能陷入某些異常情況。

ax = az.plot_trace(fit, var_names=["alpha","beta","sigma"])

Thick-haired chains ! I omitted alpha and beta[1:2] due to image size.

· 我們參數的后可信區間-森林圖。 這些應該作為我們模型參數(和超參數！)的比例和位置的直觀指南。 例如，σ(此處為PC框架[4]下的復雜度參數)不應低于0，而如果?個預測變量的可信區間不包含0，或者如果0，則我們可以了解"統計意義" 幾乎不在里面

axes = az.plot_forest( post_data, kind="forestplot", var_names= ["beta","sigma"], combined=True, ridgeplot_overlap=1.5, colors="blue", figsize=(9, 4),)

Looks good. alpha here can be seen as "collector" representing our reference category (I used reference encoding). We can certainly normalize our variables for friendlier scaling — I encourage you to play with this.

預測

貝葉斯推斷具有預測能力，我們可以通過從預測后驗分布中采樣來產生它們：

這些(以及整個貝葉斯框架)的優點在于，我們可以使用(預測的)可信區間來了解估計的方差/波動性。 畢竟，如果預測模型的預測"足夠接近"目標50的1倍，那么預測模型有什么用？

讓我們看看我們的預測值如何與測試集中的觀察值進行比較：

The straight dotted line indicates perfect prediction. We're not too far from it.

然后，我們可以對這些預測進行ML標準度量(MSE)，以根據保留的測試集中的實際值評估其質量。

當我們更改一些模型輸入時，我們還可以可視化我們的預測值(以及它們相應的95％可信度區間)：

az.style.use("arviz-darkgrid")sns.relplot(x="weight", y="mpg") data=pd.DataFrame({'weight':X_test['weight'],'mpg':y_test}))az.plot_hpd(X_test['weight'], la['y_new'], color="k", plot_kwargs={"ls": "--"})sns.relplot(x="acceleration", y="mpg", data=pd.DataFrame({'acceleration':X_test['acceleration'],'mpg':y_test}))az.plot_hpd(X_test['acceleration'], la['y_new'], color="k", plot_kwargs={"ls": "--"})

在下面，我使用統計模型將貝葉斯模型的性能與普通最大似然線性回歸模型的性能進行比較：

They perform painfully close (for this particular seed).

最佳實踐：為了更好地了解模型性能，您應該通過對K≥30的火車測試拆分運行實現，在測試MSE上引導95％置信區間(CLT)。

注意事項

· 請注意，回歸分析是許多現代數據科學問題家族的一個非常容易獲得的起點。 在學習完本教程后，我希望您開始考慮它的兩種風格：ML(基于損失)和經典(基于模型/似然性)。

· 貝葉斯推理并不難融入您的DS實踐中。可以采用基于損失的方法作為貝葉斯方法，但這是我稍后將要討論的主題。

· 只要您知道自己在做什么就很有用！ 事先指定是困難的。 但是，當您的數據量有限(昂貴的數據收集/標記，罕見事件等)或您確切知道要測試或避免的內容時，此功能特別有用。

· 在ML任務(預測，分類等)中，尤其是隨著數據集規模的增長，貝葉斯框架很少比其(頻繁)ML對手表現更好。 處理大型數據集(大數據)時，您的先驗數據很容易被淹沒。

· 正確指定的貝葉斯模型是一個生成模型，因此您應該能夠輕松生成數據，以檢查模型與相關數據集/數據生成過程的一致性。

· 在模型構建和檢查過程之前，整個過程以及之后，EDA和繪圖都是至關重要的。 [5]是一篇關于貝葉斯工作流程中可視化的重要性的出色論文。

進一步指示

我鼓勵您不僅從貝葉斯的角度而且從機器學習的角度來批判性地考慮改善此預測的方法。

數據集是否足夠大或足夠豐富，可以從嚴格的ML方法中受益？ 有什么方法可以改善這種貝葉斯模型？ 在哪種情況下，該模型肯定會勝過MLE模型？ 如果觀測值在某種程度上是相關的(聚類，自相關，空間相關等)，該怎么辦？ 當可解釋性是建模優先級(監管方面)時該怎么辦？

如果您想知道將貝葉斯方法擴展到深度學習的方法，還可以研究變分推理[6]方法，作為純貝葉斯方法的可擴展替代方法。 這是一個"簡單"的概述，并且是技術評論。

參考文獻

[1] B. Carpenter等。 Stan：一種概率編程語言(2017)。 統計軟件雜志76(1)。 DOI 10.18637 / jss.v076.i01。

[2] Betancourt。 哈密爾頓蒙特卡洛概念介紹(2017)。 arXiv：1701.02434

[3] J·韋克菲爾德。 貝葉斯和頻率回歸方法(2013)。 統計中的Springer系列。 紐約施普林格。 doi：10.1007 / 978–1–4419–0925–1。

[4] D. Simpson等。 懲罰模型組件的復雜性：構建先驗的有原則的實用方法(2017)。 在：統計科學32.1，第1至28頁。 doi：10.1214 / 16-STS576。

[5] J. Gabry等。 貝葉斯工作流中的可視化(2019)。 J. R. Stat。 Soc。 A，182：389-402。 doi：10.1111 / rssa.12378

[6] D.M. Blei等。 變異推理：《統計學家評論》(2016年)。 美國統計協會雜志，第一卷。 第112章 518，2017 DOI：10.1080 / 01621459.2017.1285773

·

(本文翻譯自Sergio E. Betancourt的文章《Painless Introduction to Applied Bayesian Inference using (Py)Stan》，參考：https://towardsdatascience.com/painless-introduction-to-applied-bayesian-inference-using-py-stan-36b503a4cd80)

總結

以上是生活随笔為你收集整理的ad20如何导入库_一文看懂如何使用（Py）Stan进行贝叶斯推理的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。