PyMC3（现在简称为PyMC）是一个贝叶斯建模包，它使数据科学家能够轻松地执行贝叶斯推理。PyMC3使用马尔可夫链蒙特卡洛(MCMC)方法计算后验分布。这种方法相当复杂。我们不会在这里详细描述原理。这里我们只解释一些简单的概念。为什么要使用MCMC？这是为了避免贝叶斯定理中计算归一化常数的棘手问题：其中P(H|D)是后验，P(H)是先验，P(D|H)是似然，P(D)是归一化常数，定义为：对于许多问题，这个积分要么没有封闭形式的解，要么无法计算。这就是为什么开发像MCMC这样的方法来解决这个问题并允许我们使用贝叶斯方法的原因。此外，还有一种叫做ConjugatePriors的方法也可以解决这个问题，但是它的可扩展性不如MCMC。如果您想了解更多有关共轭先验的信息，我们将在后面的其他文章中进行解释。在这篇文章中，我们将展示如何使用PyMC3包实现贝叶斯线性回归，并快速介绍它与普通线性回归的区别。贝叶斯回归与常客回归常客回归和贝叶斯回归方法之间的主要区别在于它们处理参数的方式。在频率统计中，线性回归模型的参数是固定的，而在贝叶斯统计中，它们是随机变量。频率论者使用最大似然估计（MLE）的方法来推导线性回归模型的值。MLE的结果是每个参数的固定值。在贝叶斯世界中，参数是具有一定概率的值分布，用更多的数据更新这个分布可以让我们更加确定参数可以取的值。此过程称为贝叶斯更新。通过上面的简要介绍，我们已经看到了贝叶斯回归和频率回归之间的主要区别。让我们言归正传，首先使用PyMC3导入包：.apiassmfif要安装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，数据=数据）。fit()results.paramsinter=results.params['Intercept']slope=results.params['x']x_vals=np.arange(min(x),max(x),0.1)ols_line=inter+斜率*x_valsfig,ax=plt.subplots(figsize=(9,5))ax.scatter(data['x'],data['y'])ax.plot(x_vals,ols_line,label='OLSFit',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个样本进行建模：#Startourmodelwithpm.Model()asmodel_100:grad=pm.Uniform("grad",lower=results.params['x']*0.5,upper=results.params['x']*1.5)inter=pm.Uniform("inter",lower=results.params['截距']*0.5,upper=results.params['截距']*1.5)sigma=pm.Uniform("sigma",lower=results.resid.std()*0.5,\upper=results.resid.std()*1.5)#线性回归线mean=inter+grad*data['x']#描述我们的条件输出的分布y=pm.Normal('y',mu=mean,sd=sigma,observed=data['y'])#使用pymc3对100个样本运行采样trace_100=pm.sample(100,return_inferencedata=True)此代码将运行MCMC采样器以计算每个参数的后验值，绘制每个参数的后验分布：withmodel_100:az.plot_posterior(trace_100,var_names=['grad','inter','sigma'],textsize=18,point_estimate='mean',rope_color='black')可以看出这些后验分布的均值和OLS估计一样，但是对于贝叶斯回归不是参数可以采用的唯一值。这里有很多值，是贝叶斯线性回归的主要核心之一。HDI代表高密度区间，它描述了我们在参数估计中的确定性。该模拟仅使用了数据中的100个样本。与其他方法一样，数据越多，贝叶斯方法就越确定。现在我们用10000个样本再次运行这个过程：withpm.Model()asmodel_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['截距']*0.5,upper=results.params['截距']*1.5)sigma=pm.Uniform("西格玛",lower=results.resid.std()*0.5,upper=results.resid.std()*1.5)#线性回归线mean=inter+grad*data['x']#描述我们的条件输出y的分布=pm.Normal('y',mu=mean,sd=sigma,observed=data['y'])#使用pymc3对10,000个样本运行采样trace_10_000=pm.sample(10_000,return_inferencedata=True)请看参数的后验分布：withmodel_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作者：EgorHowell
