複製鏈接
請複製以下鏈接發送給好友

貝葉斯線性迴歸

鎖定
貝葉斯線性迴歸(Bayesian linear regression)是使用統計學中貝葉斯推斷(Bayesian inference)方法求解的線性迴歸(linear regression)模型 [1-2] 
貝葉斯線性迴歸將線性模型的參數視為隨機變量(random variable),並通過模型參數(權重係數)的先驗(prior)計算其後驗(posterior)。貝葉斯線性迴歸可以使用數值方法求解,在一定條件下,也可得到解析型式的後驗或其有關統計量 [3] 
貝葉斯線性迴歸具有貝葉斯統計模型的基本性質,可以求解權重係數的概率密度函數,進行在線學習以及基於貝葉斯因子(Bayes factor)的模型假設檢驗 [1]  [3] 
中文名
貝葉斯線性迴歸
外文名
Bayesian linear regression
類    型
迴歸算法
提出者
Dennis Lindley,Adrian Smith
提出時間
1972年
學    科
統計學

目錄

貝葉斯線性迴歸歷史

貝葉斯線性迴歸是二十世紀60-70年代貝葉斯理論興起時得到發展的統計方法之一,其早期工作包括在迴歸模型中對權重先驗 [4-5]  和最大後驗密度(Highest Posterior Density, HPD)的研究 [6]  、在貝葉斯視角下發展的隨機效應模型(random effect mode) [7]  以及貝葉斯統計中可交換性(exchangeability)概念的引入 [8] 
正式提出貝葉斯迴歸方法的工作來自英國統計學家Dennis Lindley和Adrian Smith,在其1972年發表的兩篇論文中,Lindley和Smith對貝葉斯線性迴歸進行了系統論述並通過數值試驗與其它線性迴歸方法進行了比較,為貝葉斯線性迴歸的應用奠定了基礎 [9]  [10] 

貝葉斯線性迴歸理論與算法

貝葉斯線性迴歸模型

給定相互獨立的N組學習樣本
,貝葉斯線性迴歸使用如下的多元線性迴歸模型
[2] 
這裏
權重係數
殘差。由於學習樣本間相互獨立,因此
獨立同分布(independent and identically distributed, iid)。貝葉斯線性迴歸假設殘差服從正態分佈 [2]  ,其方差服從逆Gamma分佈(Inverse-Gamma distribution) [11] 
正態分佈的均值
和逆Gamma分佈的係數
要求預先指定。通常設定均值
,對應白噪聲(white noise)殘差,因此貝葉斯線性迴歸的模型本身至少包含2個超參數(hyperparameter)。以上的貝葉斯線性迴歸也可推廣至廣義線性模型(Generalized Linear Model, GLM),得到貝葉斯廣義線性模型(Bayesian GLM) [9] 

貝葉斯線性迴歸求解

根據線性模型的定義,權重係數
與觀測數據
相互獨立,也與殘差的方差
相互獨立,由貝葉斯定理(Bayes' theorem)可推出,貝葉斯線性迴歸中權重係數的後驗有如下表示 [2] 
式中
稱為似然(likelihood),由線性迴歸模型完全決定,以模型殘差服從iid的0均值正態分佈為例,這裏似然也服從iid的正態分佈 [2] 
式中的
的邊緣似然(marginal likelihood),在貝葉斯推斷中也被稱為模型證據(model evidence),僅與觀測數據有關,與權重係數相互獨立 [2]  [3] 
求解
式要求預先給定權重係數的先驗
,即一個連續概率分佈(continuous probability distribution),通常的選擇為0均值的正態分佈 [2] 
一維貝葉斯線性迴歸:(a,c,d)先驗、似然和後驗,(b)迴歸結果 一維貝葉斯線性迴歸:(a,c,d)先驗、似然和後驗,(b)迴歸結果
式中
為預先給定的超參數。和其它貝葉斯推斷一樣,根據似然和先驗的類型,可用於求解貝葉斯線性迴歸的方法有三類,即極大後驗估計(Maximum A Posteriori estimation, MAP)、共軛先驗(conjugate prior)求解和數值方法,前兩者要求似然或先驗滿足特定條件,數值方法沒有特定要求,可以通過迭代逼近任意形式的後驗分佈。
在權重係數沒有合理先驗假設的問題中,貝葉斯線性迴歸可使用無信息先驗(uninformative prior),即一個均勻分佈(uniform distribution),此時權重係數按均等的機會取任意值 [12] 
極大後驗估計(Maximum A Posteriori estimation, MAP)
在貝葉斯線性迴歸中,MAP可以被視為一種特殊的貝葉斯估計(Bayesian estimator) [13]  ,其求解步驟與極大似然估計(maximum likelihood estimation)類似 [13]  。對給定的先驗,MAP將
式轉化為求解
使後驗概率最大的優化問題,並求得後驗的眾數(mode)。由於正態分佈的眾數即是均值,因此MAP通常被應用於正態先驗 [2] 
這裏以的0均值正態先驗為例介紹MAP的求解步驟,首先給定權重係數
的0均值正態分佈先驗:
。由於邊緣似然與
相互獨立,此時求解後驗概率的極大值等價於求解似然和先驗乘積的極大值 [3] 
對上述極值問題取自然對數並考慮正態分佈的解析型式,可有如下推導 [3] 
由於式中各項的係數均為負數,因此該極大值問題可轉化為僅與
有關的極小值問題,並可通過線性代數得到
的解 [14] 
式中
為模型殘差和權重係數方差的比值,由超參數直接計算,
為單位矩陣。
除正態先驗外,MAP也使用拉普拉斯分佈(Laplace distribution)作為權重係數的先驗:
式中
為超參數。取均值
,在經過與正態先驗類似的推導後,拉普拉斯先驗的MAP可轉化為如下優化問題 [14] 
該優化問題對應一個泛定方程,在
足夠大時,可以得到稀疏解(sparse solution)。
MAP是單點估計,不支持在線學習(online learning),也無法提供置信區間,但能夠以很小的計算量求解貝葉斯線性迴歸的權重係數,且可用於最常見的正態先驗情形 [3]  。使用正態先驗和MAP求解的貝葉斯線性迴歸等價於嶺迴歸(ridge regression),優化目標中的
被稱為L2正則化項 [14]  ;使用拉普拉斯分佈先驗的情形對應線性模型的LASSO(Least Absolute Shrinkage and Selection Operator),
被稱為L1正則化項 [14] 
共軛先驗求解
由於貝葉斯線性迴歸的似然是正態分佈,因此在權重係數的先驗存在共軛分佈時可利用共軛性(conjugacy)求解後驗 [1]  。這裏以正態先驗為例介紹其求解步驟。
首先引入權重係數的0均值正態先驗:
,隨後由
式可知,後驗正比於似然和先驗的乘積 [2] 
在正態似然下,方差已知的正態先驗的共軛分佈是正態分佈,因此將上式按正態分佈的解析型式進行整理有如下結果 [2] 
式中
定義與先前相同。以上式為基礎,可以得到權重係數的均值和置信區間,完成對貝葉斯線性迴歸的求解。
除正態先驗外,共軛先驗求解也適用於對數正態分佈(log-normal distribution)、Beta分佈(Beta distribution)、Gamma分佈(Gamma distribution)等的先驗。
數值方法
一般地,貝葉斯推斷的數值方法都適用於貝葉斯線性迴歸,其中最常見的是馬爾可夫鏈蒙特卡羅(Markov Chain Monte Carlo, MCMC) [15]  。這裏以MCMC中的吉布斯採樣(Gibbs sampling)算法為例介紹。
給定均值為
,方差為
的正態先驗
和權重係數
,吉布斯採樣是一個迭代算法,每個迭代都依次採樣所有的權重係數 [11]  [16] 
1. 隨機初始化權重係數
和殘差的方差
2. 採樣
3. 使用
採樣
4. 採樣
5. 重複1-4至迭代完成/分佈收斂
迭代步驟中的
被稱為條件採樣分佈(conditional sampling distribution)。條件採樣的推導較為繁瑣,這裏給出一維情形下截距
斜率
和殘差方差
的條件採樣分佈 [16] 
在Python 3的Numpy模塊下,上述迭代步驟可有如下編程實現:
import numpy as np
import 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算法 [17]  和數據增強算法(data augmentation algorithm) [18]  也可用於貝葉斯線性迴歸的MCMC計算。

貝葉斯線性迴歸預測

由MAP求解的貝葉斯線性迴歸可直接使用權重係數對預測數據進行計算。MAP是單點估計,因此預測結果不提供後驗分佈。對共軛先驗或數值方法求解的貝葉斯線性迴歸,可通過邊緣化(marginalized out)模型權重,即按其後驗積分得到預測結果 [2] 
式中
為預測數據,
為預測結果。對0均值的正態先驗,由其共軛性可知,預測結果的後驗也為正態分佈 [2] 
上式右側的正態分佈可以提供預測結果的置信區間,在本質上是線性模型使用所有潛在權重係數計算所得結果的概率密度函數。
模型驗證
邊緣似然描述了似然和先驗的組合在多大程度上解釋了觀測數據,因此邊緣似然可以用於計算貝葉斯因子,與其它模型進行比較並驗證模型和先驗的合理性 [1]  。由全概率公式(law of total probability)可知,邊緣似然有如下積分形式 [2] 
由上式計算的邊緣似然擁有描述模型複雜度的全部信息,包括先驗假設的類型,權重係數的數量等,因此可以與任意複雜度的其它模型計算貝葉斯因子進行比較。
由MCMC計算的貝葉斯線性迴歸也可以使用數值方法進行交叉驗證(cross-validation),具體方法包括重採樣(re-sampling MCMC) [19]  和使用EC(Expectation Consistent)近似的留一法(Leave One Out, LOO)交叉驗證 [20] 
在線學習
由共軛先驗和數值方法求解的貝葉斯線性迴歸可以進行在線學習,即使用實時數據對權重係數進行更新。在線學習的具體方法依模型本身而定,其設計思路是將先錢求解得到的後驗作為新的先驗,並帶入數據得到後驗 [3] 

貝葉斯線性迴歸性質

穩健性:由求解部分的推導可知,若貝葉斯線性迴歸使用正態先驗,則其MAP的估計結果等價於嶺迴歸,而使用拉普拉斯先驗的情形對應線性模型的LASSO,因此貝葉斯線性迴歸與使用正則化(regularization)的迴歸分析一樣平衡了模型的經驗風險和結構風險。特別地,使用拉普拉斯先驗的貝葉斯線性迴歸由於可以得到稀疏解,因此具有一定的變量篩選(variable selection)能力 [21] 
作為貝葉斯推斷所具有的性質:由於貝葉斯線性迴歸是貝葉斯推斷在線性迴歸問題中的應用,因此其具有貝葉斯方法的一般性質,包括對先驗進行實時更新、將觀測數據視為定點因而不需要漸進假設、服從似然定理(likelihood principle)、估計結果包含置信區間 [1]  。基於最小二乘法的線性迴歸(Ordinary Linear Regression, OLR)通常僅在觀測數據顯著地多於權重係數維數的時候才會有好的效果,而貝葉斯線性迴歸沒有此類限制,且在權重係數維數過高是,可以根據後驗對模型進行降維(dimensionality reduction) [22] 
與高斯過程迴歸(Gaussian Process Regression, GPR)的比較:貝葉斯線性迴歸是GPR在權重空間(weight-space)下的特例。在貝葉斯線性迴歸中引入映射函數
可得到GPR的預測形式:
,而當映射函數為等值函數(identity function)時,GPR退化為貝葉斯線性迴歸 [2] 

貝葉斯線性迴歸應用

除一般意義上線性迴歸模型的應用外,貝葉斯線性模型可被用於觀測數據較少但要求提供後驗分佈的問題,例如對物理常數的精確估計 [23]  。有研究利用貝葉斯線性迴歸的性質進行變量篩選 [21]  和降維 [22]  。此外,貝葉斯線性迴歸是在統計方法中使用貝葉斯推斷的簡單實現之一,因此常作為貝葉斯理論或數值計算教學的重要例子 [16] 
參考資料
  • 1.    Wakefield, J..Bayesian and frequentist regression methods:Springer Science & Business Media,2013:Chapter 3
  • 2.    Rasmussen, C.E. and Williams, C.K.I..Gaussian processes in machine learning:MIT Press,2006:pp.7-11
  • 3.    Polykovskiy, D. and Novikov, A., Bayesian Methods for Machine Learning  .Coursera and National Research University Higher School of Economics.2017[引用日期2018-11-23]
  • 4.    Fisk, P.R., 1967. Models of the second kind in regression analysis. Journal of the Royal Statistical Society. Series B (Methodological), pp.266-281.
  • 5.    Nelder, J.A., 1968. Regression, model-building and invariance. Journal of the Royal Statistical Society. Series A (General), pp.303-329.
  • 6.    Box, G.E. and Tiao, G.C., 1965. Multiparameter problems from a Bayesian point of view. The Annals of Mathematical Statistics, 36(5), pp.1468-1482.
  • 7.    Box, G.E. and Tiao, G.C., 1968. Bayesian estimation of means for the random effect model. Journal of the American Statistical Association, 63(321), pp.174-181.
  • 8.    De Finetti, B., 1964. Foresight: Its logical laws, its subjective sources. Studies in subjective probability, 1964, pp.94-158.
  • 9.    Smith, A.F., 1973. A general Bayesian linear model. Journal of the Royal Statistical Society. Series B (Methodological), pp.67-75.
  • 10.    Lindley, D.V. and Smith, A.F., 1972. Bayes estimates for the linear model. Journal of the Royal Statistical Society. Series B (Methodological), pp.1-41.
  • 11.    Wand, M., Markov Chain Monte Carlo for Bayesian Linear Regression  .School of Mathematical and Physical Sciences, University of Technology Sydney[引用日期2018-11-23]
  • 12.    Stauffer, H. B..Contemporary Bayesian and Frequentist Statistical Research Methods for Natural Resource Scientists :John Wiley and Sons, Inc.,2008:Chapter 2
  • 13.    Bassett, R. and Deride, J., 2016. Maximum a posteriori estimators as a limit of Bayes estimators. Mathematical Programming, pp.1-16.
  • 14.    Rai, P., Learning via Probabilistic Modeling, Probabilistic Linear Regression  .Computer Science and Engineering, Indian Institute of Technology (IIT) Kanpur.2016-8-12[引用日期2018-11-23]
  • 15.    Evans, S., 2012. Bayesian regression analysis. Electronic Theses and Dissertations, University of Louisville (Paper 421)
  • 16.    Neal, R. M., CSC 2541: Bayesian Methods for Machine Learning  .Department of Computer Science, University of Toronto.2011[引用日期2018-11-22]
  • 17.    Kalinova, V., Lecture notes on Regression: Markov Chain Monte Carlo (MCMC)   .Machine Learning course: the elegant way to extract information from data, Max Planck Institute for Radioastronomy.2017-2-13[引用日期2018-11-23]
  • 18.    Choi, H.M. and Hobert, J.P., 2013. Analysis of MCMC algorithms for Bayesian linear regression with Laplace errors. Journal of Multivariate Analysis, 117, pp.32-40.
  • 19.    Bhattacharya, S. and Haslett, J., 2007. Importance re-sampling MCMC for cross-validation in inverse problems. Bayesian Analysis, 2(2), pp.385-407.
  • 20.    Kabashima, Y., Obuchi, T. and Uemura, M., 2016. Approximate cross-validation formula for Bayesian linear regression. In Communication, Control, and Computing (Allerton), 2016 54th Annual Allerton Conference on (pp. 596-600). IEEE.
  • 21.    Mitchell, T.J. and Beauchamp, J.J., 1988. Bayesian variable selection in linear regression. Journal of the American Statistical Association, 83(404), pp.1023-1032.
  • 22.    Bernardo, J.M., Bayarri, M.J., Berger, J.O., Dawid, A.P., Heckerman, D., Smith, A.F.M. and West, M., 2003. Bayesian factor regression models in the “large p, small n” paradigm. Bayesian statistics, 7, pp.733-742.
  • 23.    Wübbeler, G., Bodnar, O. and Elster, C., 2017. Robust Bayesian linear regression with application to an analysis of the CODATA values for the Planck constant. Metrologia, 55(1), p.20.
展開全部 收起