Python用PyMC3实现贝叶斯线性回归模型

原文链接：http://tecdat.cn/?p=5263

在本文中，我们将在贝叶斯框架中引入回归建模，并使用PyMC3 MCMC库进行推理。

我们将首先回顾经典或频率论者的多重线性回归方法。然后我们将讨论贝叶斯如何考虑线性回归。

用PyMC3进行贝叶斯线性回归

在本节中，我们将对统计实例进行一种历史悠久的方法，即模拟一些我们知道的属性的数据，然后拟合一个模型来恢复这些原始属性。

什么是广义线性模型？

在我们开始讨论贝叶斯线性回归之前，我想简要地概述广义线性模型（GLM）的概念，因为我们将使用它们来在PyMC3中制定我们的模型。

广义线性模型是将普通线性回归扩展到更一般形式的回归的灵活机制，包括逻辑回归（分类）和泊松回归（用于计数数据）以及线性回归本身。

GLM允许具有除正态分布以外的误差分布的响应变量（参见频率分区中的上述）。

用PyMC3模拟数据并拟合模型

在我们使用PyMC3来指定和采样贝叶斯模型之前，我们需要模拟一些噪声线性数据。

if __name__ == "__main__":
    # 参数真值
    beta_0 = 1.0  # Intercept
    beta_1 = 2.0  # Slope
 
    #模拟
    N = 100
    eps_sigma_sq = 0.5
 
    # 模拟
    df = simulate_linear_data(N, beta_0, beta_1, eps_sigma_sq)
 
    # 可视化
    sns.lmplot(x="x", y="y", data=df, size=10)
    plt.xlim(0.0, 1.0)

输出如下图所示：

通过Numpy，pandas和seaborn模拟噪声线性数据

现在我们已经进行了模拟，我们想要对数据拟合贝叶斯线性回归。这是glm模块进来的地方。它使用与R指定模型类似的模型规范语法。

然后我们将找到MCMC采样器的最大后验概率（MAP）估计值。最后，我们将使用No-U-Turn Sampler（NUTS）来进行实际推理，然后绘制模型的曲线，将前500个样本丢弃为“burn in”

traceplot如下图所示：

使用PyMC3将贝叶斯GLM线性回归模型拟合到模拟数据

我们可以使用glm库调用的方法绘制这些线plot_posterior_predictive。该方法采用trace对象和plot（samples）的行数。

首先我们使用seaborn lmplot方法，这次fit_reg参数设置False为停止绘制频数回归曲线。然后我们绘制100个采样的后验预测回归线。最后，我们绘制使用原始的“真实”回归线和β1=2的参数。下面的代码片段产生了这样的情节：β0=1β0=1β1=2β1=2

我们可以在下图中看到回归线的抽样范围：