這是為了避開(kāi)貝葉斯定理中計(jì)算歸一化常數(shù)的棘手問(wèn)題：

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

對(duì)于許多問(wèn)題，這個(gè)積分要么沒(méi)有封閉形式的解，要么無(wú)法計(jì)算。所以才有MCMC等方法被開(kāi)發(fā)出來(lái)解決這個(gè)問(wèn)題，并允許我們使用貝葉斯方法。
此外還有一種叫做共軛先驗(yàn)(Conjugate Priors)的方法也能解決這個(gè)問(wèn)題，但它的可延展性不如MCMC。如果你想了解更多關(guān)于共軛先驗(yàn)的知識(shí)，我們?cè)诤竺嫫渌恼逻M(jìn)行講解。
在這篇文章中，我們將介紹如何使用PyMC3包實(shí)現(xiàn)貝葉斯線性回歸，并快速介紹它與普通線性回歸的區(qū)別。

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

首先導(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，請(qǐng)查看他們網(wǎng)站上的安裝說(shuō)明。現(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()

我們先使用頻率論的常見(jiàn)回歸，使用普通最小二乘(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)行對(duì)比。
要使用PyMC3，我們必須初始化一個(gè)模型，選擇先驗(yàn)并告訴模型后驗(yàn)分布應(yīng)該是什么，我們使用100個(gè)樣本來(lái)進(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)

該代碼將運(yùn)行MCMC采樣器來(lái)計(jì)算每個(gè)參數(shù)的后驗(yàn)值，繪制每個(gè)參數(shù)的后驗(yàn)分布：

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

可以看到這些后驗(yàn)分布的平均值與OLS估計(jì)相同，但對(duì)于貝葉斯回歸來(lái)說(shuō)并不是參數(shù)可以采用的唯一值。這里有很多值，這是貝葉斯線性回歸的主要核心之一。HDI代表高密度區(qū)間(High Density Interval)，它描述了我們?cè)趨?shù)估計(jì)中的確定性。
這個(gè)模擬只使用了數(shù)據(jù)中的100個(gè)樣本。和其他方法一樣，數(shù)據(jù)越多，貝葉斯方法就越確定。現(xiàn)在我們使用10000個(gè)樣本，再次運(yùn)行這個(gè)過(guò)程：

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ù)的后驗(yàn)分布：

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

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

來(lái)源：DeepHub IMBA