在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<lower=0> N;      int<lower=0,upper=1> y[N];    }    parameters {      real<lower=0,upper=1> 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<lower=1> N; // number of training samples    int<lower=0> K; // number of predictors - 1 (intercept)    matrix[N, K] x; // matrix of predictors    vector[N] y_obs; // observed/training mpg        int<lower=1> N_new;    matrix[N_new, K] x_new;}parameters {    real alpha;    vector[K] beta;    real<lower=0> 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（）？ 数学统计中有一个众所周知的结果，该结果表明，长度为<∞且具有单位对角协方差矩阵且向量均值μ<∞的多重正态随机变量（向量W）等效于以下独立独立正态随机变量w的向量 N（μ，1）分布。

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》，参考：)