在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）筆記本並開始我們的工作。 首先，讓我們使用以下內容導入庫：

<code>import pystan 

import pickle
import numpy as np
import arviz as az
import pandas as pd
import seaborn as sns
import statsmodels.api as statmod
import matplotlib.pyplot as plt
from IPython.display import Image
from IPython.core.display import HTML/<code>

投幣推理

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

我們選擇了以下模型：

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

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

<code>sns.distplot(np.random.beta(5,5, size=10000),kde=False)/<code>

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

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

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

· 模型對應於我們的先驗（β）和可能性（bernoulli）的定義。

<code># bernoulli model
model_code = """
    data {
      int<lower> N;
      int<lower> y[N];
    }
    parameters {
      real<lower> 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 arrays
print(fit.stansummary())/<lower>/<lower>/<lower>/<code>

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

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

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

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

以下是上述模型的輸出：

<code>

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

應用迴歸：每加侖汽車和英里數（MPG）

Image by Susan Sewert from Pixabay

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

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

<code>

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

現在，對於貝葉斯公式：

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

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

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

· 我的可能性是多少？

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

· 我的參數（和超參數）是什麼？應該選擇哪種先驗條件？

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

EDA

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

首先加載數據集：

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

讓我們檢查一下目標變量和預測變量之間的關係（如果有）：

<code>

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

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

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

<code>

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.

訓練/擬合

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

<code>from numpy import random
from sklearn import preprocessing, metrics, linear_model
from sklearn.model_selection import train_test_split
from sklearn.metrics import mean_squared_error
random.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 category
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.15, random_state=0)
X_train.head()/<code>

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

<code># Succinct matrix notation
cars_code = """
data {
    int<lower> N; // number of training samples
    int<lower> K; // number of predictors - 1 (intercept)
    matrix[N, K] x; // matrix of predictors
    vector[N] y_obs; // observed/training mpg
    
    int<lower> N_new;
    matrix[N_new, K] x_new;
}
parameters {
    real alpha;
    vector[K] beta;
    real<lower> 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
}
"""/<lower>/<lower>/<lower>/<lower>/<code>

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

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

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

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

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

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

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

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

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

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

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

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

<code>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 model
la = fit.extract(permuted=True)
print(fit.stansummary())/<code>

這是我打印的輸出（出於本教程的目的，我擺脫了That和n_eff）：

<code>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 
...
...  /<code>

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

診斷

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

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

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

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

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

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

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倍，那麼預測模型有什麼用？

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

<code>

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

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

<code>

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

<code>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": "--"})/<code>

<code>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": "--"})/<code>

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

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)

閱讀更多 聞數起舞 的文章

關鍵字: 筆記本電腦 社交網絡 推理