貝葉斯線性回歸

貝葉斯線性回歸（Bayesian linear regression）是使用統計學中貝葉斯推斷（Bayesian inference）方法求解的線性回歸（linear regression）模型。

貝葉斯線性回歸將線性模型的參數視為隨機變數（random variable），並通過模型參數（權重係數）的先驗（prior）計算其後驗（posterior）。貝葉斯線性回歸可以使用數值方法求解，在一定條件下，也可得到解析型式的後驗或其有關統計量。

貝葉斯線性回歸具有貝葉斯統計模型的基本性質，可以求解權重係數的機率密度函式，進行線上學習以及基於貝葉斯因子（Bayes factor）的模型假設檢驗。

基本介紹

• 中文名：貝葉斯線性回歸
• 外文名：Bayesian linear regression
• 類型：回歸算法
• 提出者：Dennis Lindley，Adrian Smith
• 提出時間：1972年
• 學科：統計學

歷史

貝葉斯線性回歸是二十世紀60-70年代貝葉斯理論興起時得到發展的統計方法之一，其早期工作包括在回歸模型中對權重先驗和最大後驗密度（Highest Posterior Density, HPD）的研究、在貝葉斯視角下發展的隨機效應模型（random effect mode）以及貝葉斯統計中可交換性（exchangeability）概念的引入。

正式提出貝葉斯回歸方法的工作來自英國統計學家Dennis Lindley和Adrian Smith，在其1972年發表的兩篇論文中，Lindley和Smith對貝葉斯線性回歸進行了系統論述並通過數值試驗與其它線性回歸方法進行了比較，為貝葉斯線性回歸的套用奠定了基礎。

理論與算法

模型

給定相互獨立的N組學習樣本 ， ，貝葉斯線性回歸使用如下的多元線性回歸模型 ：

這裡 為權重係數， 為殘差。由於學習樣本間相互獨立，因此 為獨立同分布（independent and identically distributed, iid）。貝葉斯線性回歸假設殘差服從常態分配，其方差服從逆Gamma分布（Inverse-Gamma distribution）：

常態分配的均值 和逆Gamma分布的係數 要求預先指定。通常設定均值 ，對應白噪聲（white noise）殘差，因此貝葉斯線性回歸的模型本身至少包含2個超參數（hyperparameter）。以上的貝葉斯線性回歸也可推廣至廣義線性模型（Generalized Linear Model, GLM），得到貝葉斯廣義線性模型（Bayesian GLM）。

求解

根據線性模型的定義，權重係數 與觀測數據 相互獨立，也與殘差的方差 相互獨立，由貝葉斯定理（Bayes' theorem）可推出，貝葉斯線性回歸中權重係數的後驗有如下表示：

式中 稱為似然（likelihood），由線性回歸模型完全決定，以模型殘差服從iid的0均值常態分配為例，這裡似然也服從iid的常態分配：

式中的 是 的邊緣似然（marginal likelihood），在貝葉斯推斷中也被稱為模型證據（model evidence），僅與觀測數據有關，與權重係數相互獨立。

求解 式要求預先給定權重係數的先驗 ，即一個連續機率分布（continuous probability distribution），通常的選擇為0均值的常態分配：

式中 為預先給定的超參數。和其它貝葉斯推斷一樣，根據似然和先驗的類型，可用於求解貝葉斯線性回歸的方法有三類，即極大後驗估計（Maximum A Posteriori estimation, MAP）、共軛先驗（conjugate prior）求解和數值方法，前兩者要求似然或先驗滿足特定條件，數值方法沒有特定要求，可以通過疊代逼近任意形式的後驗分布。

一維貝葉斯線性回歸：（a,c,d）先驗、似然和後驗，（b）回歸結果

在權重係數沒有合理先驗假設的問題中，貝葉斯線性回歸可使用無信息先驗（uninformative prior），即一個均勻分布（uniform distribution），此時權重係數按均等的機會取任意值。

極大後驗估計（Maximum A Posteriori estimation, MAP）

在貝葉斯線性回歸中，MAP可以被視為一種特殊的貝葉斯估計（Bayesian estimator），其求解步驟與極大似然估計（maximum likelihood estimation）類似。對給定的先驗，MAP將 式轉化為求解 使後驗機率最大的最佳化問題，並求得後驗的眾數（mode）。由於常態分配的眾數即是均值，因此MAP通常被套用於正態先驗。

這裡以的0均值正態先驗為例介紹MAP的求解步驟，首先給定權重係數 的0均值常態分配先驗： 。由於邊緣似然與 相互獨立，此時求解後驗機率的極大值等價於求解似然和先驗乘積的極大值：

對上述極值問題取自然對數並考慮常態分配的解析型式，可有如下推導：

由於式中各項的係數均為負數，因此該極大值問題可轉化為僅與 有關的極小值問題，並可通過線性代數得到 的解：

式中 為模型殘差和權重係數方差的比值，由超參數直接計算， 為單位矩陣。

除正態先驗外，MAP也使用拉普拉斯分布（Laplace distribution）作為權重係數的先驗：

式中 為超參數。取均值 ，在經過與正態先驗類似的推導後，拉普拉斯先驗的MAP可轉化為如下最佳化問題：

該最佳化問題對應一個泛定方程，在 足夠大時，可以得到稀疏解（sparse solution）。

MAP是單點估計，不支持線上學習（online learning），也無法提供置信區間，但能夠以很小的計算量求解貝葉斯線性回歸的權重係數，且可用於最常見的正態先驗情形。使用正態先驗和MAP求解的貝葉斯線性回歸等價於嶺回歸（ridge regression），最佳化目標中的 被稱為L₂正則化項；使用拉普拉斯分布先驗的情形對應線性模型的LASSO（Least Absolute Shrinkage and Selection Operator）， 被稱為L₁正則化項。

共軛先驗求解

由於貝葉斯線性回歸的似然是常態分配，因此在權重係數的先驗存在共軛分布時可利用共軛性（conjugacy）求解後驗。這裡以正態先驗為例介紹其求解步驟。

首先引入權重係數的0均值正態先驗： ，隨後由 式可知，後驗正比於似然和先驗的乘積：

在正態似然下，方差已知的正態先驗的共軛分布是常態分配，因此將上式按常態分配的解析型式進行整理有如下結果：

式中 定義與先前相同。以上式為基礎，可以得到權重係數的均值和置信區間，完成對貝葉斯線性回歸的求解。

除正態先驗外，共軛先驗求解也適用於對數常態分配（log-normal distribution）、Beta分布（Beta distribution）、Gamma分布（Gamma distribution）等的先驗。

數值方法

一般地，貝葉斯推斷的數值方法都適用於貝葉斯線性回歸，其中最常見的是馬爾可夫鏈蒙特卡羅（Markov Chain Monte Carlo, MCMC）。這裡以MCMC中的吉布斯採樣（Gibbs sampling）算法為例介紹。

給定均值為 ，方差為 的正態先驗 和權重係數 ，吉布斯採樣是一個疊代算法，每個疊代都依次採樣所有的權重係數：

1. 隨機初始化權重係數 和殘差的方差

2. 採樣 ：

3. 使用 採樣 ：

4. 採樣 ：

5. 重複1-4至疊代完成/分布收斂

疊代步驟中的 和 被稱為條件採樣分布（conditional sampling distribution）。條件採樣的推導較為繁瑣，這裡給出一維情形下截距 、斜率 和殘差方差 的條件採樣分布：

在Python 3的Numpy模組下，上述疊代步驟可有如下編程實現：

import numpy as npimport matplotlib.pyplot as plt# 構建測試數據## 自變數：[0,5]區間均勻分布x = np.random.uniform(low=0, high=5, size=5)## 因變數：y=2x-1+eps, eps=norm(0, 0.5)y = np.random.normal(2*x-1, 0.5)def calculate_csd(x, y, w1, w2, s, N, hyper_param):    '''    條件採樣分布函式: (w1, w2, s)        x , y: 輸入數據        w1, w2: 線性回歸模型的截距和斜率        s: 殘差方差        N: 樣本量    '''    # 超參數    mu1, s1 = hyper_param["mu1"], hyper_param['s1']    mu2, s2 = hyper_param["mu2"], hyper_param['s2']    a, b, = hyper_param["a"], hyper_param["b"]    # 計算方差    var1 = (s1*s)/(s+s1*N)    var2 = (s*s2)/(s+s2*np.sum(x*x))    # 計算均值    mean1 = (s*mu1+s1*np.sum(y-w2*x))/(s+s1*N)    mean2 = (s*mu2+s2*np.sum((y-w1)*x))/(s+s2*np.sum(x*x))    # 計算(a, b)    eps = y-w1-w2*x    a = a+N/2; b = b+np.sum(eps*eps)/2    # 使用numpy.random採樣（inv-gamma=1/gamma）    # 注意numpy.random中normal和gamma的定義方式：    ##     np.random.normal(mean, std)    ##     inv-Gamma = 1/np.random.gamma(a, scale), scale=1/b    return np.random.normal(mean1, np.sqrt(var1)), \           np.random.normal(mean2, np.sqrt(var2)), \           1/np.random.gamma(a, 1/b)def Gibbs_sampling(x, y, iters, init, hyper_param):    '''    吉布斯採樣算法主程式        iter: 疊代次數        init: 初始值        hyper_param: 超參數    '''    N = len(y) # 樣本量    # 採樣初始值    w1, w2, s = init["w1"], init["w2"], init["s"]    # 使用數組記錄疊代值    history = np.zeros([iters, 3])*np.nan    # for循環疊代    for i in range(iters):        w1, w2, s = calculate_csd(x, y, w1, w2, s, N, hyper_param)        history[i, :] = np.array([w1, w2, s])            return history# 開始採樣iters = 10000 # MCMC疊代1e4次init = {'w1': 0, 'w2': 0, 's': 0}hyper_param = {'mu1': 0, 's1': 1, 'mu2': 0, 's2': 1, 'a': 2, 'b': 1}history = Gibbs_sampling(x, y, iters, init, hyper_param)burnt = history[500:] # 截取收斂後的馬爾可夫鏈ax1 = plt.subplot(2, 1, 1); ax2 = plt.subplot(2, 1, 2)ax1.hist(burnt[:, 0], bins=100)ax1.text(2, 0.025*iters, '$\mathrm{N(w_1|\mu,s)}$', fontsize=14)ax2.hist(burnt[:, 1], bins=100)ax2.text(.5, 0.025*iters, '$\mathrm{N(w_2|\mu,s)}$', fontsize=14)

除吉布斯採樣外，Metropolis-Hastings算法和數據增強算法（data augmentation algorithm）也可用於貝葉斯線性回歸的MCMC計算。

預測

由MAP求解的貝葉斯線性回歸可直接使用權重係數對預測數據進行計算。MAP是單點估計，因此預測結果不提供後驗分布。對共軛先驗或數值方法求解的貝葉斯線性回歸，可通過邊緣化（marginalized out）模型權重，即按其後驗積分得到預測結果：

式中 為預測數據， 為預測結果。對0均值的正態先驗，由其共軛性可知，預測結果的後驗也為常態分配：

上式右側的常態分配可以提供預測結果的置信區間，在本質上是線性模型使用所有潛在權重係數計算所得結果的機率密度函式。

模型驗證

邊緣似然描述了似然和先驗的組合在多大程度上解釋了觀測數據，因此邊緣似然可以用於計算貝葉斯因子，與其它模型進行比較並驗證模型和先驗的合理性。由全機率公式（law of total probability）可知，邊緣似然有如下積分形式：

由上式計算的邊緣似然擁有描述模型複雜度的全部信息，包括先驗假設的類型，權重係數的數量等，因此可以與任意複雜度的其它模型計算貝葉斯因子進行比較。

由MCMC計算的貝葉斯線性回歸也可以使用數值方法進行交叉驗證（cross-validation），具體方法包括重採樣（re-sampling MCMC）和使用EC（Expectation Consistent）近似的留一法（Leave One Out, LOO）交叉驗證。

線上學習

由共軛先驗和數值方法求解的貝葉斯線性回歸可以進行線上學習，即使用實時數據對權重係數進行更新。線上學習的具體方法依模型本身而定，其設計思路是將先錢求解得到的後驗作為新的先驗，並帶入數據得到後驗。

性質

穩健性：由求解部分的推導可知，若貝葉斯線性回歸使用正態先驗，則其MAP的估計結果等價於嶺回歸，而使用拉普拉斯先驗的情形對應線性模型的LASSO，因此貝葉斯線性回歸與使用正則化（regularization）的回歸分析一樣平衡了模型的經驗風險和結構風險。特別地，使用拉普拉斯先驗的貝葉斯線性回歸由於可以得到稀疏解，因此具有一定的變數篩選（variable selection）能力。

作為貝葉斯推斷所具有的性質：由於貝葉斯線性回歸是貝葉斯推斷線上性回歸問題中的套用，因此其具有貝葉斯方法的一般性質，包括對先驗進行實時更新、將觀測數據視為定點因而不需要漸進假設、服從似然定理（likelihood principle）、估計結果包含置信區間等。基於最小二乘法的線性回歸（Ordinary Linear Regression, OLR）通常僅在觀測數據顯著地多於權重係數維數的時候才會有好的效果，而貝葉斯線性回歸沒有此類限制，且在權重係數維數過高是，可以根據後驗對模型進行降維（dimensionality reduction）。

與高斯過程回歸（Gaussian Process Regression, GPR）的比較：貝葉斯線性回歸是GPR在權重空間（weight-space）下的特例。在貝葉斯線性回歸中引入映射函式 可得到GPR的預測形式： ，而當映射函式為等值函式（identity function）時，GPR退化為貝葉斯線性回歸。

套用

除一般意義上線性回歸模型的套用外，貝葉斯線性模型可被用於觀測數據較少但要求提供後驗分布的問題，例如對物理常數的精確估計。有研究利用貝葉斯線性回歸的性質進行變數篩選和降維。此外，貝葉斯線性回歸是在統計方法中使用貝葉斯推斷的簡單實現之一，因此常作為貝葉斯理論或數值計算教學的重要例子。