在PyStan中應(yīng)用貝葉斯回歸

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

設(shè)定

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

以下是設(shè)置Stan的詳細(xì)安裝步驟：

對(duì)于MacOS：

· 安裝miniconda / anaconda

· 安裝xcode

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

· 創(chuàng)建一個(gè)環(huán)境stan或pystan

· 輸入conda install numpy來安裝numpy或?qū)umpy替換為您需要安裝的軟件包

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

· 或者：pip install pystan

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

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

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

投幣推理

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

我們選擇了以下模型：

下面的概率質(zhì)量函數(shù)對(duì)應(yīng)于Y：

首先，通過使用NumPy的隨機(jī)功能進(jìn)行仿真，我們可以了解參數(shù)為a = b = 5的先驗(yàn)行為：

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

如概念背景中所述，此先驗(yàn)似乎是合理的，因?yàn)槠鋵?duì)稱性約為0.5，但仍可能在任一方向上產(chǎn)生偏差。 然后，我們可以使用以下語法在pystan中定義模型：

· 數(shù)據(jù)對(duì)應(yīng)于我們模型的數(shù)據(jù)。 在這種情況下，整數(shù)N對(duì)應(yīng)于拋硬幣的次數(shù)，y對(duì)應(yīng)于長(zhǎng)度為N的整數(shù)的向量，該向量將包含我的實(shí)驗(yàn)觀察結(jié)果。

· 參數(shù)對(duì)應(yīng)于我們模型的參數(shù)，在這種情況下為theta或獲得"頭部"的概率。

· 模型對(duì)應(yīng)于我們的先驗(yàn)(β)和可能性(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())

請(qǐng)注意，model.sampling的默認(rèn)參數(shù)是iter = 1000，chains = 4和warmup = 500。 我們會(huì)根據(jù)時(shí)間和可用的計(jì)算資源進(jìn)行調(diào)整。

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

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

· chains≥0對(duì)應(yīng)于我們樣本中的MCMC鏈數(shù)。

以下是上述模型的輸出：

· 我們?chǔ)鹊暮缶导s為0.36 <0.5。 盡管如此，95％的后可信區(qū)間很寬：(0.14，0.61)。 因此，我們可以認(rèn)為結(jié)果在統(tǒng)計(jì)上不是結(jié)論性的，但它指向偏見而沒有跳到0。

應(yīng)用回歸：每加侖汽車和英里數(shù)(MPG)

Image by Susan Sewert from Pixabay

讓我們構(gòu)建一個(gè)貝葉斯線性回歸模型來解釋和預(yù)測(cè)汽車數(shù)據(jù)集中的MPG。 盡管我的方法是傳統(tǒng)的基于可能性的模型，但我們可以(并且應(yīng)該！)使用原始ML訓(xùn)練檢驗(yàn)分裂來評(píng)估我們的預(yù)測(cè)質(zhì)量。

該數(shù)據(jù)集來自UCI ML存儲(chǔ)庫(kù)，包含以下信息：

我們?yōu)槭裁床粓?jiān)持我們的基礎(chǔ)知識(shí)并使用標(biāo)準(zhǔn)的線性回歸？ [3]回顧其以下功能形式：

現(xiàn)在，對(duì)于貝葉斯公式：

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

您可以在此處詢問有關(guān)Σ的信息。 這是一個(gè)N(訓(xùn)練觀察數(shù))x P(系數(shù)/特征數(shù))矩陣。 在標(biāo)準(zhǔn)線性回歸情況下，要求此Σ是對(duì)角矩陣，且沿對(duì)角線有σ(觀測(cè)值之間的獨(dú)立性)。

在執(zhí)行EDA時(shí)，您應(yīng)該始終考慮以下幾點(diǎn)：

· 我的可能性是多少？

· 我的模特應(yīng)該是什么？ 互動(dòng)與沒有互動(dòng)？ 多項(xiàng)式？

· 我的參數(shù)(和超參數(shù))是什么？應(yīng)該選擇哪種先驗(yàn)條件？

· 我應(yīng)該考慮任何聚類，時(shí)間或空間依賴性嗎？

EDA

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

首先加載數(shù)據(jù)集：

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

讓我們檢查一下目標(biāo)變量和預(yù)測(cè)變量之間的關(guān)系(如果有)：

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

現(xiàn)在，讓我們檢查一下數(shù)字變量和mpg之間的關(guān)系：

· 我更喜歡形象化這些關(guān)系而不是計(jì)算相關(guān)性，因?yàn)樗o了我關(guān)于它們之間關(guān)系的視覺和數(shù)學(xué)意義，超越了簡(jiǎn)單的標(biāo)量。

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.

訓(xùn)練/擬合

讓我們準(zhǔn)備數(shù)據(jù)集以進(jìn)行擬合和測(cè)試：

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()

現(xiàn)在進(jìn)入我們的模型規(guī)范，它與上面的拋硬幣問題沒有實(shí)質(zhì)性的區(qū)別，除了我使用矩陣符號(hào)來盡可能簡(jiǎn)化模型表達(dá)式的事實(shí)之外：

# 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}"""

· 數(shù)據(jù)與我們模型的數(shù)據(jù)相對(duì)應(yīng)，如上面的拋硬幣示例所示。

· 參數(shù)對(duì)應(yīng)于我們模型的參數(shù)。

· 此處經(jīng)過轉(zhuǎn)換的參數(shù)使我可以將theta定義為模型在訓(xùn)練集上的擬合值

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

在對(duì)初始數(shù)據(jù)進(jìn)行檢查之后，選擇了以上先驗(yàn)條件：

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

· 那么α和β呢？ 在正常情況下，最簡(jiǎn)單的方法是為這些參數(shù)選擇常規(guī)先驗(yàn)。

· 為什么要對(duì)β使用multi_normal()？ 數(shù)學(xué)統(tǒng)計(jì)中有一個(gè)眾所周知的結(jié)果，該結(jié)果表明，長(zhǎng)度為

stan的有趣之處在于，我們可以要求在測(cè)試集上進(jìn)行預(yù)測(cè)而無需重新擬合-而是可以在第一個(gè)調(diào)用中從stan請(qǐng)求它們。

· 我們通過在上面的單元格中定義t_new來實(shí)現(xiàn)這一點(diǎn)。

· theta是我們對(duì)訓(xùn)練集的合適預(yù)測(cè)。

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

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()是一個(gè)像表一樣排列的字符串，它為我提供了擬合期間估計(jì)的每個(gè)參數(shù)的后驗(yàn)均值，標(biāo)準(zhǔn)差和幾個(gè)百分點(diǎn)。 雖然alpha對(duì)應(yīng)于截距，但我們有8個(gè)beta回歸變量，復(fù)雜度或觀察誤差sigma，訓(xùn)練集theta的擬合值以及測(cè)試集y_new的預(yù)測(cè)值。

診斷

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

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

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

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

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

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.

預(yù)測(cè)

貝葉斯推斷具有預(yù)測(cè)能力，我們可以通過從預(yù)測(cè)后驗(yàn)分布中采樣來產(chǎn)生它們：

這些(以及整個(gè)貝葉斯框架)的優(yōu)點(diǎn)在于，我們可以使用(預(yù)測(cè)的)可信區(qū)間來了解估計(jì)的方差/波動(dòng)性。 畢竟，如果預(yù)測(cè)模型的預(yù)測(cè)"足夠接近"目標(biāo)50的1倍，那么預(yù)測(cè)模型有什么用？

讓我們看看我們的預(yù)測(cè)值如何與測(cè)試集中的觀察值進(jìn)行比較：

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

然后，我們可以對(duì)這些預(yù)測(cè)進(jìn)行ML標(biāo)準(zhǔn)度量(MSE)，以根據(jù)保留的測(cè)試集中的實(shí)際值評(píng)估其質(zhì)量。

當(dāng)我們更改一些模型輸入時(shí)，我們還可以可視化我們的預(yù)測(cè)值(以及它們相應(yīng)的95％可信度區(qū)間)：

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

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

They perform painfully close (for this particular seed).

最佳實(shí)踐：為了更好地了解模型性能，您應(yīng)該通過對(duì)K≥30的火車測(cè)試拆分運(yùn)行實(shí)現(xiàn)，在測(cè)試MSE上引導(dǎo)95％置信區(qū)間(CLT)。

注意事項(xiàng)

· 請(qǐng)注意，回歸分析是許多現(xiàn)代數(shù)據(jù)科學(xué)問題家族的一個(gè)非常容易獲得的起點(diǎn)。 在學(xué)習(xí)完本教程后，我希望您開始考慮它的兩種風(fēng)格：ML(基于損失)和經(jīng)典(基于模型/似然性)。

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

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

· 在ML任務(wù)(預(yù)測(cè)，分類等)中，尤其是隨著數(shù)據(jù)集規(guī)模的增長(zhǎng)，貝葉斯框架很少比其(頻繁)ML對(duì)手表現(xiàn)更好。 處理大型數(shù)據(jù)集(大數(shù)據(jù))時(shí)，您的先驗(yàn)數(shù)據(jù)很容易被淹沒。

· 正確指定的貝葉斯模型是一個(gè)生成模型，因此您應(yīng)該能夠輕松生成數(shù)據(jù)，以檢查模型與相關(guān)數(shù)據(jù)集/數(shù)據(jù)生成過程的一致性。

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

進(jìn)一步指示

我鼓勵(lì)您不僅從貝葉斯的角度而且從機(jī)器學(xué)習(xí)的角度來批判性地考慮改善此預(yù)測(cè)的方法。

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

如果您想知道將貝葉斯方法擴(kuò)展到深度學(xué)習(xí)的方法，還可以研究變分推理[6]方法，作為純貝葉斯方法的可擴(kuò)展替代方法。 這是一個(gè)"簡(jiǎn)單"的概述，并且是技術(shù)評(píng)論。

參考文獻(xiàn)

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

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

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

[4] D. Simpson等。 懲罰模型組件的復(fù)雜性：構(gòu)建先驗(yàn)的有原則的實(shí)用方法(2017)。 在：統(tǒng)計(jì)科學(xué)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等。 變異推理：《統(tǒng)計(jì)學(xué)家評(píng)論》(2016年)。 美國(guó)統(tǒng)計(jì)協(xié)會(huì)雜志，第一卷。 第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)

總結(jié)

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

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯(cuò)，歡迎將生活随笔推薦給好友。