关于机器学习:贝叶斯回归使用-PyMC3-实现贝叶斯回归

PyMC3(当初简称为 PyMC)是一个贝叶斯建模包，它使数据科学家可能轻松地进行贝叶斯推断。

PyMC3 采纳马尔可夫链蒙特卡罗 (MCMC) 办法计算后验散布。这个办法相当简单，原理方面咱们这里不做详细描述，这里只阐明一些简略的概念，为什么应用 MCMC 呢？

这是为了避开贝叶斯定理中计算归一化常数的辣手问题:

其中 P(H | D)为后验，P(H)为先验，P(D | H)为似然，P(D)为归一化常数，定义为:

对于许多问题，这个积分要么没有关闭模式的解，要么无奈计算。所以才有 MCMC 等办法被开发进去解决这个问题，并容许咱们应用贝叶斯办法。

此外还有一种叫做共轭先验 (Conjugate Priors) 的办法也能解决这个问题，但它的可延展性不如 MCMC。如果你想理解更多对于共轭先验的常识，咱们在前面其余文章进行解说。

在这篇文章中，咱们将介绍如何应用 PyMC3 包实现贝叶斯线性回归，并疾速介绍它与一般线性回归的区别。

贝叶斯 vs 频率回归

频率主义和贝叶斯回归办法之间的要害区别在于他们如何解决参数。在频率统计中，线性回归模型的参数是固定的，而在贝叶斯统计中，它们是随机变量。

频率主义者应用极大似然预计 (MLE) 的办法来推导线性回归模型的值。MLE 的后果是每个参数的一个固定值。

在贝叶斯世界中，参数是具备肯定概率的值散布，应用更多的数据更新这个散布，这样咱们就能够更加确定参数能够取的值。这个过程被称为贝叶斯更新

有了下面的简略介绍，咱们曾经晓得了贝叶斯和频率回归之间的次要区别。上面开始正题

应用 PyMC3

首先导入包：

 
 import pymc3 as pm
 import arviz as az
 
 import pandas as pd
 import numpy as np
 import matplotlib.pyplot as plt
 from sklearn import datasets
 from scipy.stats import norm
 import statsmodels.formula.api as smf

如果想装置 PyMC3 和 ArviZ，请查看他们网站上的装置阐明。当初咱们应用 sklearn 中的 make_regression 函数生成一些数据:

 x, y = datasets.make_regression(n_samples=10_000,
                                 n_features=1,
                                 noise=10,
                                 bias=5)
 
 data = pd.DataFrame(list(zip(x.flatten(), y)),columns =['x', 'y'])
 
 fig, ax = plt.subplots(figsize=(9,5))
 ax.scatter(data['x'], data['y'])
 ax.ticklabel_format(style='plain')
 plt.xlabel('x',fontsize=18)
 plt.ylabel('y',fontsize=18)
 plt.xticks(fontsize=18)
 plt.yticks(fontsize=18)
 plt.show()

咱们先应用频率论的常见回归，应用一般最小二乘 (OLS) 办法绘制频率线性回归线:

 formula = 'y ~ x'
 results = smf.ols(formula, data=data).fit()
 results.params
 
 inter = results.params['Intercept']
 slope = results.params['x']
 x_vals = np.arange(min(x), max(x), 0.1)
 ols_line = inter + slope * x_vals
 
 fig, ax = plt.subplots(figsize=(9,5))
 ax.scatter(data['x'], data['y'])
 ax.plot(x_vals, ols_line,label='OLS Fit', color='red')
 ax.ticklabel_format(style='plain')
 plt.xlabel('x',fontsize=18)
 plt.ylabel('y',fontsize=18)
 plt.xticks(fontsize=18)
 plt.yticks(fontsize=18)
 plt.legend(fontsize=16)
 plt.show()

下面的后果咱们作为基线模型与咱们前面的贝叶斯回归进行比照

要应用 PyMC3，咱们必须初始化一个模型，抉择先验并通知模型后验散布应该是什么，咱们应用 100 个样本来进行建模，:

 # Start our model
 with pm.Model() as model_100:
     grad = pm.Uniform("grad",
                       lower=results.params['x']*0.5,
                       upper=results.params['x']*1.5)
     
     inter = pm.Uniform("inter",
                        lower=results.params['Intercept']*0.5,
                        upper=results.params['Intercept']*1.5)
     
     sigma = pm.Uniform("sigma",
                        lower=results.resid.std()*0.5,\
                        upper=results.resid.std()*1.5)
     
     # Linear regression line
     mean = inter + grad*data['x']
     
     # Describe the distribution of our conditional output
     y = pm.Normal('y', mu = mean, sd = sigma, observed = data['y'])
     # Run the sampling using pymc3 for 100 samples
     trace_100 = pm.sample(100,return_inferencedata=True)

该代码将运行 MCMC 采样器来计算每个参数的后验值，绘制每个参数的后验散布:

 with model_100:
     az.plot_posterior(trace_100,
                       var_names=['grad', 'inter', 'sigma'],
                       textsize=18,
                       point_estimate='mean',
                       rope_color='black')

能够看到这些后验散布的平均值与 OLS 预计雷同，但对于贝叶斯回归来说并不是参数能够采纳的惟一值。这里有很多值，这是贝叶斯线性回归的次要外围之一。HDI 代表高密度区间(High Density Interval)，它形容了咱们在参数估计中的确定性。

这个模仿只应用了数据中的 100 个样本。和其余办法一样，数据越多，贝叶斯办法就越确定。当初咱们应用 10000 个样本，再次运行这个过程：

 with pm.Model() as model_10_100: 
     grad = pm.Uniform("grad",
                       lower=results.params['x']*0.5,
                       upper=results.params['x']*1.5)
     
     inter = pm.Uniform("inter",
                        lower=results.params['Intercept']*0.5,
                        upper=results.params['Intercept']*1.5)
     
     sigma = pm.Uniform("sigma",
                        lower=results.resid.std()*0.5,
                        upper=results.resid.std()*1.5)
     
     # Linear regression line
     mean = inter + grad*data['x']
     
     # Describe the distribution of our conditional output
     y = pm.Normal('y', mu = mean, sd = sigma, observed = data['y'])
     # Run the sampling using pymc3 for 10,000 samples
     trace_10_000 = pm.sample(10_000,return_inferencedata=True)

看看参数的后验散布:

 with model_10_100:
     az.plot_posterior(trace_10_000,
                       var_names=['grad', 'inter', 'sigma'],
                       textsize=18,
                       point_estimate='mean',
                       rope_color='black')

能够看到平均值的预测并没有扭转，然而随着对参数的散布更加确定，总体上的散布变得更加平滑和紧凑。

总结

在本文中，咱们介绍贝叶斯统计的次要原理，并解释了它与频率统计相比如何采纳不同的办法进行线性回归。而后，咱们学习了如何应用 PyMC3 包执行贝叶斯回归的根本示例。

本文的代码：

https://avoid.overfit.cn/post/07ce671c5022406a8d299bfa196af871

作者：Egor Howell