這是為了避開貝葉斯定理中計算歸一化常數(shù)的棘手問題：

其中P(H | D)為后驗，P(H)為先驗，P(D | H)為似然，P(D)為歸一化常數(shù)，定義為：

對于許多問題，這個積分要么沒有封閉形式的解，要么無法計算。所以才有MCMC等方法被開發(fā)出來解決這個問題，并允許我們使用貝葉斯方法。
此外還有一種叫做共軛先驗(Conjugate Priors)的方法也能解決這個問題，但它的可延展性不如MCMC。如果你想了解更多關(guān)于共軛先驗的知識，我們在后面其他文章進(jìn)行講解。
在這篇文章中，我們將介紹如何使用PyMC3包實現(xiàn)貝葉斯線性回歸，并快速介紹它與普通線性回歸的區(qū)別。

頻率主義和貝葉斯回歸方法之間的關(guān)鍵區(qū)別在于他們?nèi)绾翁幚韰?shù)。在頻率統(tǒng)計中，線性回歸模型的參數(shù)是固定的，而在貝葉斯統(tǒng)計中，它們是隨機變量。
頻率主義者使用極大似然估計(MLE)的方法來推導(dǎo)線性回歸模型的值。MLE的結(jié)果是每個參數(shù)的一個固定值。
在貝葉斯世界中，參數(shù)是具有一定概率的值分布，使用更多的數(shù)據(jù)更新這個分布，這樣我們就可以更加確定參數(shù)可以取的值。這個過程被稱為貝葉斯更新。
有了上面的簡單介紹，我們已經(jīng)知道了貝葉斯和頻率回歸之間的主要區(qū)別。下面開始正題：

首先導(dǎo)入包：

 









 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，請查看他們網(wǎng)站上的安裝說明?，F(xiàn)在我們使用sklearn中的make_regression函數(shù)生成一些數(shù)據(jù)：

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

上面的結(jié)果我們作為基線模型與我們后面的貝葉斯回歸進(jìn)行對比。
要使用PyMC3，我們必須初始化一個模型，選擇先驗并告訴模型后驗分布應(yīng)該是什么，我們使用100個樣本來進(jìn)行建模：

 # 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采樣器來計算每個參數(shù)的后驗值，繪制每個參數(shù)的后驗分布：

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

可以看到這些后驗分布的平均值與OLS估計相同，但對于貝葉斯回歸來說并不是參數(shù)可以采用的唯一值。這里有很多值，這是貝葉斯線性回歸的主要核心之一。HDI代表高密度區(qū)間(High Density Interval)，它描述了我們在參數(shù)估計中的確定性。
這個模擬只使用了數(shù)據(jù)中的100個樣本。和其他方法一樣，數(shù)據(jù)越多，貝葉斯方法就越確定。現(xiàn)在我們使用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)

看看參數(shù)的后驗分布：

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

可以看到平均值的預(yù)測并沒有改變，但是隨著對參數(shù)的分布更加確定，總體上的分布變得更加平滑和緊湊。
總結(jié)
在本文中，我們介紹貝葉斯統(tǒng)計的主要原理，并解釋了它與頻率統(tǒng)計相比如何采用不同的方法進(jìn)行線性回歸。然后，我們學(xué)習(xí)了如何使用PyMC3包執(zhí)行貝葉斯回歸的基本示例。
本文的代碼：https://github.com/egorhowell/Medium-Articles/blob/main/Statistics/pymc3_tutorial.ipynb
作者：Egor Howell

來源：DeepHub IMBA