貝葉斯線性回歸（Bayesian Linear Regression）

2016年06月21日 09:50:40 Duanxx 閱讀數 54254更多

分類專欄： 監督學習

 

版權聲明：本文為博主原創文章，遵循 CC 4.0 BY-SA 版權協議，轉載請附上原文出處鏈接和本聲明。

本文鏈接： https://blog.csdn.net/daunxx/article/details/51725086

貝葉斯線性回歸（Bayesian Linear Regression）

標簽（空格分隔）： 監督學習

---

@ author : [e-mail] 
@ time : 2015-06-19

原文地址

 

• 貝葉斯線性回歸Bayesian Linear Regression
  • 原文地址
  • 關於參數估計
  • 極大似然估計
    • 漸進無偏
    • 漸進一致
  • 最大后驗估計
  • 貝葉斯估計
    • 貝葉斯估計核心問題
    • 貝葉斯估計第一個重要元素
    • 貝葉斯估計第二個重要元素
    • 貝葉斯估計的增量學習
  • 貝葉斯線性回歸
    • 貝葉斯線性回歸的學習過程
    • 貝葉斯回歸的優缺點
    • 貝葉斯脊回歸Bayesian Ridge Regression

 

本文的研究順序是：

極大似然估計最大后驗估計貝葉斯估計貝葉斯線性回歸

關於參數估計

在很多的機器學習或數據挖掘的問題中，我們所面對的只有數據，但數據中潛在的概率密度函數是不知道的，其概率密度分布需要我們從數據中估計出來。想要確定數據對應的概率密度分布，就需要確定兩個東西：概率密度函數的形式 和 概率密度函數的參數。

有時可能知道的是概率密度函數的形式(高斯、瑞利等等)，但是不知道具體的參數，例如均值或者方差；還有的時候可能不知道概率密度的類型，但是知道一些估計的參數，比如均值和方差。

關於上面提到的需要確定的兩個東西：概率密度函數的形式和參數，至少在機器學習的教課書上，我所看到的情況都是：給了一堆數據，然后假設其概率密度函數的形式為 高斯分布 ，或者是混合高斯分布，那么，剩下的事情就是對高斯分布的參數，μμ 和 σ2σ2 進行估計。所以，參數估計，便成了極其最重要的問題。

其實，常用的參數估計方法有：極大似然估計、最大后驗估計、貝葉斯估計、最大熵估計、混合模型估計。他們之間是有遞進關系的，想要理解后一個參數估計方法，最好對前一個參數估計有足夠的理解。

要想清晰的說明貝葉斯線性回歸，或者叫做貝葉斯參數估計，就必須對極大似然估計、最大后驗估計做詳細的說明，他們之間是有遞進的關系的。

極大似然估計

在之前《多項式回歸》的文章中，用最后一小節是線性回歸的概率解釋，其中就說明了以平方誤差維損失函數的最小二乘法和極大似然估計的等價性，在這個基礎上，本文更為詳細的討論極大似然估計。

這里先以一個分類問題來說明一般參數估計面對的數據形式。考慮一個MM類的問題，特征向量服從p(x|ωi),i=1,2...,Mp(x|ωi),i=1,2...,M 分布。這是現實情況中最常見的一種數據存在形式，數據集合XX是由MM個類別的數據子集Xm,m=1,2...,MXm,m=1,2...,M 組成的，第mm類別的數據子集XmXm 對應的概率密度函數是p(x|ωm)p(x|ωm)。

前面已經介紹過了，想要確定數據的概率分布，需要知道概率密度函數的 形式 和 參數，這里首先做一個基本假設：概率分布的形式已知，比如假設每個類別的數據都滿足高斯分布，那么，似然函數就可以以參數  θi  θi 的形式表示，如果是高斯分布，則參數為μiμi和σ2iσi2，即θi=(μi,σ2i)θi=(μi,σi2)。

為了強調概率分布p(x|ωi)p(x|ωi)和 θiθi 有關，將對應的概率密度函數記為p(x|ωi;θi)p(x|ωi;θi)，這種記法屬於頻率概率學派的記法。這里的極大似然估計對應於一個類條件概率密度函數。

在概率論中一直有兩大學派，分別是頻率學派和貝葉斯學派。簡單點說，頻率學派認為，概率是頻率的極限，比如投硬幣，當實驗次數足夠大時，正面朝上的頻率可以認為是這枚硬幣正面朝上的概率，這個是頻率學派。但是，如果要預測一些未發生過的事情，比如，北極的冰山在2050年完全融化的概率，由於這個事情完全沒有發生過，所以無法用頻率來代替概率表示，只能研究過去幾十年，北極冰山融化的速率，並將其作為先驗條件，來預測北極的冰山在2050年完全融化的概率，這就是概率的貝葉斯學派。上面的問題，如果用貝葉斯學派的記法的話，是：p(x|ωi,θi)p(x|ωi,θi)。這兩個學派適用的情況不太一樣，但是，在我目前所用到的概率論的知識中，貌似這兩個學派並沒有什么太大的區別，只是記法略有不同，稍微注意下即可。

從上面的描述中可以知道，利用每一個類XiXi中已知的特征向量集合，可以估計出其對應的參數θiθi。進一步假設每一類中的數據不影響其他類別數據的參數估計，那么上面的MM個類別的參數估計就可以用下面這個統一的模型，獨立的解決：

設x1,x2,...,xNx1,x2,...,xN 是從概率密度函數p(x;θ)p(x;θ)中隨機抽取的樣本，那么就可以得到聯合概率密度函數 p(X;θ)p(X;θ)， 其中X={x1,x2,...,xN}X={x1,x2,...,xN} 是樣本集合。假設不同的樣本之間具有統計獨立性，那么：

 

 

p(X;θ)≡p(x1,x2,...,xN;θ)=∏k=1Np(xk;θ)p(X;θ)≡p(x1,x2,...,xN;θ)=∏k=1Np(xk;θ)

 

注意：這里的p(xk;θ)p(xk;θ) 本來的寫法是 p(x|ωi;θi)p(x|ωi;θi) , 是一個類條件概率密度函數，只是因為這里是一個統一的模型，所以可以將wiwi 省略。

需要重申一下，想要得到上面這個公式，是做了幾個基本的假設的，第一：假設MM個類別的數據子集的概率密度函數形式一樣，只是參數的取值不同；第二：假設類別ii中的數據和類別jj中的數據是相互獨立抽樣的，即類別jj的參數僅僅根據類別jj的數據就可以估計出來，類別ii的數據並不能影響類別jj的參數估計，反之亦然；第三：每個類別內的樣本之間具有統計獨立性，即每個類別內的樣本之間是獨立同分布 (iidiid) 的。

此時，就可以使用最大似然估計(Maximum Likelihood,ML)來估計參數θθ了：

 

 

θ^ML=argmaxθ∏k=1Np(xk;θ)θ^ML=argmaxθ∏k=1Np(xk;θ)

 

為了得到最大值，θ^MLθ^ML必須滿足的必要條件是,似然函數對θθ的梯度必須為0，即：

 

 

∂∏Nk=1p(xk;θ)∂θ=0∂∏k=1Np(xk;θ)∂θ=0

 

一般我們取其對數形式：

 

 

L(θ)≡ln∏k=1Np(xk;θ)L(θ)≡ln∏k=1Np(xk;θ)

 

 

 

∂L(θ)∂θ=∑k=1N∂lnp(xk;θ)∂θ=∑k=1N1p(xk;θ)∂p(xk;θ)∂θ=0∂L(θ)∂θ=∑k=1N∂lnp(xk;θ)∂θ=∑k=1N1p(xk;θ)∂p(xk;θ)∂θ=0

 

需要注意：極大似然估計對應於似然函數的峰值

極大似然估計有兩個非常重要的性質：漸進無偏 和 漸進一致性，有了這兩個性質，使得極大似然估計的成為了非常簡單而且實用的參數估計方法。這里假設θ0θ0是密度函數p(x;θ)p(x;θ)中未知參數的准確值。

漸進無偏

極大似然估計是漸進無偏的，即：

 

 

limN→∞E[θ^ML]=θ0limN→∞E[θ^ML]=θ0

 

也就是說，這里認為估計值 θ^MLθ^ML 本身是一個隨機變量（因為不同的樣本集合XX會得到不同的 θ^MLθ^ML），那么其均值就是未知參數的真實值，這就是漸進無偏。

漸進一致

極大似然估計是漸進一致的，即：

 

 

limN→∞prob{∥θ^ML−θ0∥⩽ϵ}=1limN→∞prob{‖θ^ML−θ0‖⩽ϵ}=1

 

這個公式還可以表示為：

 

 

limN→∞E∥θ^ML−θ0∥2=0limN→∞E‖θ^ML−θ0‖2=0

 

對於一個估計器而言，一致性是非常重要的，因為存在滿足無偏性，但是不滿足一致性的情況，比如，θ^MLθ^ML在 θ0θ0周圍震盪。如果不滿足一致性，那么就會出現很大的方差。

注意：以上兩個性質，都是在漸進的前提下（N→∞N→∞）才能討論的，即只有NN足夠大時，上面兩個性質才能成立

最大后驗估計

在最大似然估計（MAP）中，將θθ看做是未知的參數，說的通俗一點，最大似然估計是θθ的函數，其求解過程就是找到使得最大似然函數最大的那個參數θθ。

從最大后驗估計開始，將參數θθ 看成一個隨機變量，並在已知樣本集{x1,x2,...,xN}{x1,x2,...,xN} 的條件下，估計參數θθ。

這里一定要注意，在最大似然估計中，參數θθ是一個定值，只是這個值未知，最大似然函數是θθ的函數，這里θθ是沒有概率意義的。但是，在最大后驗估計中，θθ是有概率意義的，θθ有自己的分布，而這個分布函數，需要通過已有的樣本集合XX得到，即最大后驗估計需要計算的是 p(θ|X)p(θ|X)

根據貝葉斯理論：

 

 

p(θ|X)=p(θ)p(X|θ)p(X)p(θ|X)=p(θ)p(X|θ)p(X)

 

這就是參數θθ關於已有數據集合XX的后驗概率，要使得這個后驗概率最大，和極大似然估計一樣，這里需要對后驗概率函數求導。由於分子中的p(X)p(X)相對於θθ是獨立的，隨意可以直接忽略掉p(X)p(X)。

 

 

θ^MAP=argmaxθp(θ|X)=argmaxθp(θ)p(X|θ)θ^MAP=argmaxθp(θ|X)=argmaxθp(θ)p(X|θ)

 

為了得到參數θθ，和ML一樣，需要對p(θ|X)p(θ|X)求梯度，並使其等於0：

 

 

p(θ|X)∂θ=p(θ)p(X|θ）∂θ=0p(θ|X)∂θ=p(θ)p(X|θ）∂θ=0

 

注意：這里p(X|θ)p(X|θ)和極大似然估計中的似然函數p(X;θ)p(X;θ)是一樣的，只是記法不一樣。MAP和ML的區別是：MAP是在ML的基礎上加上了p(θ)p(θ)

這里需要說明，雖然從公式上來看 MAP=ML∗p(θ)MAP=ML∗p(θ)，但是這兩種算法有本質的區別，ML將θθ視為一個確定未知的值，而MAP則將θθ視為一個隨機變量。

在MAP中，p(θ)p(θ)稱為θθ的先驗，假設其服從均勻分布，即對於所有θθ取值，p(θ)p(θ)都是同一個常量，則MAP和ML會得到相同的結果。當然了，如果p(θ)p(θ)的方差非常的小，也就是說，p(θ)p(θ)是近似均勻分布的話，MAP和ML的結果自然也會非常的相似。

貝葉斯估計

注意：以下所有的概率分布表述方式均為貝葉斯學派的表述方式。

貝葉斯估計核心問題

為了防止標號混淆，這里定義已有的樣本集合為DD，而不是之前的XX。樣本集合DD中的樣本都是從一個 固定但是未知 的概率密度函數p(x)p(x)中獨立抽取出來的，要求根據這些樣本估計xx的概率分布，記為p(x|D)p(x|D)，並且使得p(x|D)p(x|D)盡量的接近p(x)p(x)，這就是貝葉斯估計的核心問題。

貝葉斯估計第一個重要元素

雖然p(x)p(x)是未知的，但是前面提到過，一個密度分布的兩個要素為：形式和參數，我們可以假設p(x)p(x)的形式已知，但是參數θθ的取值未知。這里就有了貝葉斯估計的第一個重要元素 p(x|θ)p(x|θ)，這是一個條件概率密度函數，准確的說，是一個類條件概率密度函數（具體原因參見本文前面關於極大似然估計的說明）。強調一下：p(x|θ)p(x|θ)的形式是已知的，只是參數θθ的取值未知。由於這里的xx可以看成一個測試樣本，所以這個條件密度函數，從本質上講，是θθ 在點 xx 處的似然估計。

貝葉斯估計第二個重要元素

由於參數θθ的取值未知，且，我們將θθ看成一個隨機變量，那么，在觀察到具體的訓練樣本之前，關於θθ的全部知識，可以用一個先驗概率密度函數p(θ)p(θ)來表示，對於訓練樣本的觀察，使得我們能夠把這個先驗概率密度轉化成為后驗概率密度函數p(θ|D)p(θ|D)，根據后驗概率密度相關的論述知道，我們希望p(θ|D)p(θ|D)在θθ的真實值附近有非常顯著的尖峰。這里的這個后驗概率密度，就是貝葉斯估計的第二個主要元素。

現在，將貝葉斯估計核心問題p(x|D)p(x|D)，和貝葉斯估計的兩個重要元素：p(x|θ)p(x|θ)、p(θ|D)p(θ|D)聯系起來：

 

 

p(x|D)=∫p(x,θ|D)dθ=∫p(x|θ,D)p(θ|D)dθp(x|D)=∫p(x,θ|D)dθ=∫p(x|θ,D)p(θ|D)dθ

 

上面式子中，xx是測試樣本，DD是訓練集，xx和DD的選取是獨立進行的，因此，p(x|θ,D)p(x|θ,D)可以寫成p(x|θ)p(x|θ)。所以，貝葉斯估計的核心問題就是下面這個公式：

 

 

p(x|D)=∫p(x|θ)p(θ|D)dθp(x|D)=∫p(x|θ)p(θ|D)dθ

 

下面這句話一定要理解：這里p(x|θ)p(x|θ)是θθ關於測試樣本xx這一個點的似然估計，而p(θ|D)p(θ|D)則是θθ在已有樣本集合上的后驗概率。這就是為什么本文一開始並沒有直接講貝葉斯估計，而是先說明極大似然估計和最大后驗估計的原因。其中，后驗概率p(θ|D)p(θ|D)為：

 

 

p(θ|D)=p(D|θ)p(θ)p(D)=p(D|θ)p(θ)∫p(D|θ)p(θ)dθp(θ|D)=p(D|θ)p(θ)p(D)=p(D|θ)p(θ)∫p(D|θ)p(θ)dθ

 

 

 

p(D|θ)=∏k=1Np(xk|θ)p(D|θ)=∏k=1Np(xk|θ)

 

上面這個式子就是貝葉斯估計最核心的公式，它把類條件概率密度p(x|D)p(x|D) (這里一定要理解為什么是類條件概率密度，其實這個的准確寫法可以是p(x|Dm)p(x|Dm),或者 p(x|wm,Dm)p(x|wm,Dm)，具體原因參見本文前面關於極大似然估計的部分) 和未知參數向量θθ的后驗概率密度p(θ|D)p(θ|D)聯系在了一起。如果后延概率密度p(θ|D)p(θ|D)在某一個值θ^θ^附近形成顯著的尖峰，那么就有p(x|D)≈p(x|θ^)p(x|D)≈p(x|θ^)，就是說，可以用估計值 θ^θ^近似代替真實值所得的結果。

貝葉斯估計的增量學習

為了明確的表示樣本集合DD中有nn個樣本，這里采用記號：Dn={x1,x2,...,xn}Dn={x1,x2,...,xn}。根據前一個公式，在n>1n>1的情況下有：

 

 

p(Dn|θ)=p(xn|θ)p(Dn−1|θ)p(Dn|θ)=p(xn|θ)p(Dn−1|θ)

 

可以很容易得到：

 

 

p(θ|Dn)=p(xn|θ)p(Dn−1|θ)p(θ)∫p(xn|θ)p(Dn−1|θ)p(θ)dθ=p(xn|θ)p(θ|Dn−1)∫p(xn|θ)p(θ|Dn−1)dθp(θ|Dn)=p(xn|θ)p(Dn−1|θ)p(θ)∫p(xn|θ)p(Dn−1|θ)p(θ)dθ=p(xn|θ)p(θ|Dn−1)∫p(xn|θ)p(θ|Dn−1)dθ

 

當沒有觀測樣本時，定義p(θ|D0)=p(θ)p(θ|D0)=p(θ)，為參數θθ的初始估計。然后讓樣本集合依次進入上述公式，就可以得到一系列的概率密度函數：p(θ|D0)p(θ|D0)、p(θ|D1)p(θ|D1)、p(θ|D2)p(θ|D2)、…、p(θ|Dn)p(θ|Dn)，這一過程稱為參數估計貝葉斯遞歸法，也叫貝葉斯估計的增量學習。這是一個在線學習算法，它和隨機梯度下降法有很多相似之處。

貝葉斯線性回歸

根據之前的文章《線性回歸》、《多項式回歸》中關於極大似然估計的說明，以及本文前面關於極大似然估計的論述，可以很容易知道，如果要將極大似然估計應用到線性回歸模型中，模型的復雜度會被兩個因素所控制：基函數的數目和樣本的數目。盡管為對數極大似然估計加上一個正則項（或者是參數的先驗分布），在一定程度上可以限制模型的復雜度，防止過擬合，但基函數的選擇對模型的性能仍然起着決定性的作用。

上面說了那么大一段，就是想說明一個問題：由於極大似然估計總是會使得模型過於的復雜以至於產生過擬合的現象，所以單純的適用極大似然估計並不是特別的有效。

當然，交叉驗證是一種有效的限制模型復雜度，防止過擬合的方法，但是交叉驗證需要將數據分為訓練集合測試集，對數據樣本的浪費也是非常的嚴重的。

基於上面的討論，這里就可以引出本文的核心內容：貝葉斯線性回歸。貝葉斯線性回歸不僅可以解決極大似然估計中存在的過擬合的問題，而且，它對數據樣本的利用率是100%，僅僅使用訓練樣本就可以有效而准確的確定模型的復雜度。

這里面對的模型是線性回歸模型，其詳細的介紹可以參見前面的文章《線性回歸》,線性回歸模型是一組輸入變量xx的基函數的線性組合，在數學上其形式如下： 

 

y(x,w)=w0+∑j=1Mωjϕj(x)y(x,w)=w0+∑j=1Mωjϕj(x)

 

這里ϕj(x)ϕj(x)就是前面提到的基函數，總共的基函數的數目為MM個，如果定義ϕ0(x)=1ϕ0(x)=1的話，那個上面的式子就可以簡單的表示為： 

 

y(x,w)=∑j=0Mωjϕj(x)=wTϕ(x)y(x,w)=∑j=0Mωjϕj(x)=wTϕ(x)

 

 

 

w=(w0,w1,w2,...,wM)w=(w0,w1,w2,...,wM)

 

 

 

ϕ=(ϕ0,ϕ1,ϕ2,...,ϕM)ϕ=(ϕ0,ϕ1,ϕ2,...,ϕM)

則線性模型的概率表示如下：

 

 

 

p(t|x,w,β)=N(t|y(x,w),β−1I)p(t|x,w,β)=N(t|y(x,w),β−1I)

 

假設參數ww滿足高斯分布，這是一個先驗分布：

 

 

p(w)=N(w|0,α−1I)p(w)=N(w|0,α−1I)

 

一般來說，我們稱p(w)p(w)為共軛先驗(conjugate prior)。這里tt是xx對應的目標輸出，β−1β−1和α−1α−1分別對應於樣本集合和ww的高斯分布的方差，ww是參數，

那么，線性模型的對數后驗概率函數：

 

 

lnp(θ|D)=lnp(w|T)=−β2∑n=1N{y(xn,w)−tn}2+α2wTw+constlnp(θ|D)=lnp(w|T)=−β2∑n=1N{y(xn,w)−tn}2+α2wTw+const

 

這里TT是數據樣本的目標值向量，T={t1,t2,...,tn}T={t1,t2,...,tn}，const是和參數ww無關的量。關於這個式子的具體推導，可參見前面的文章《多項式回歸》。

貝葉斯線性回歸的學習過程

根據前面關於貝葉斯估計的增量學習可以很容易得到下面這個式子，這個就是貝葉斯學習過程：在前一個訓練集合Dn−1Dn−1的后驗概率p(θ|Dn−1)p(θ|Dn−1)上，乘以新的測試樣本點xnxn的似然估計，得到新的集合DnDn的后驗概率p(θ|Dn)p(θ|Dn)，這樣，相當於p(θ|Dn−1)p(θ|Dn−1)成為了p(θ|Dn)p(θ|Dn)的先驗概率分布：

 

 

p(θ|Dn)∝p(xn|θ)p(θ|Dn−1)p(θ|Dn)∝p(xn|θ)p(θ|Dn−1)

 

有了上面的基礎知識，這里就着重的講下面這幅圖，這個圖是從RMPL第155頁截取下來的，這幅圖清晰的描述了貝葉斯線性回歸的學習過程，下面結合這幅圖，詳細的說明一下貝葉斯學習過程。

首先，說一下這里的模型：

 

 

y(x,w)=w0+w1xy(x,w)=w0+w1x

 

第一行： 
第一行是初始狀態，此時只有關於ww的先驗信息，即：p(θ|D0)=p(θ)=N(w|0,α−1I)p(θ|D0)=p(θ)=N(w|0,α−1I)。先看中間這幅圖，這幅圖是關於ww的先驗分布，由於我們假設ww初始為高斯分布N(w|0,α−1I)N(w|0,α−1I)，故其圖形是以(0,0)(0,0)為中心的圓組成的。由於此時還沒有樣本點進入，所以沒有關於樣本的似然估計，故第一行中左邊likelihood沒有圖。第一行右邊data space的那幅圖顯示的是從第二幅圖prior/posterior中隨機抽取一些點(w0,w1)(w0,w1)，並以(w0,w1)(w0,w1)為參數，畫出來的直線，此時這些直線是隨機的。

第二行： 
此時有了第一個樣本點x1x1，那么根據x1x1就可以得到第二行中，關於x1x1的似然估計，由於y=w0+w1xy=w0+w1x，似然估計的結果其實是這個式子的對偶式，即w1=1xy−1xw0w1=1xy−1xw0。從第二行的最右邊data space中的圖中可以估計出，第一個樣本點的坐標大概為:(0.9,0.1)(0.9,0.1)，所以其第一幅圖中，似然估計的中心線的方程為:

 

 

w1=19−109w0w1=19−109w0

 

近似為左邊那幅圖的畫法。由於第二行的先驗分布是第一行的后驗分布，也就是第一行的中間那幅圖。則，第二行的后驗分布的求法就是：將第二行的第左邊那幅圖和第一行的中間那幅圖相乘，就可以得到第二行中間那幅圖。第二行最右邊那幅圖就是從第二行中間那幅圖中隨機抽取一些點(w0,w1)(w0,w1)，並以(w0,w1)(w0,w1)為參數，畫出來的直線。

第三行之后，就可以一次類推了。

上面就是貝葉斯學習過程的完整描述。

貝葉斯回歸的優缺點

優點： 
1. 貝葉斯回歸對數據有自適應能力，可以重復的利用實驗數據，並防止過擬合 
2. 貝葉斯回歸可以在估計過程中引入正則項 
缺點： 
1. 貝葉斯回歸的學習過程開銷太大

貝葉斯脊回歸（Bayesian Ridge Regression）

前面已經證明過了，如果貝葉斯線性回歸的先驗分布為

 

 

p(w)=N(w|0,α−1I)p(w)=N(w|0,α−1I)

 

那么，其最終的后驗分布公式為：

 

 

lnp(θ|D)=lnp(w|T)=−β2∑n=1N{y(xn,w)−tn}2+α2wTw+constlnp(θ|D)=lnp(w|T)=−β2∑n=1N{y(xn,w)−tn}2+α2wTw+const

 

這個相當於脊回歸，所以將這種特殊情況稱為貝葉斯脊回歸，它擁有脊回歸的所有特性，具體可以參見前面的文章《脊回歸》。

下面這份代碼提供了貝葉斯回歸的用法，以及其和最小二乘法的比較。

#!/usr/bin/python # -*- coding: utf-8 -*- """ author ： [e-mail] time : 2016-06-21-09-21 貝葉斯脊回歸 這里在一個自己生成的數據集合上測試貝葉斯脊回歸 貝葉斯脊回歸和最小二乘法(OLS)得到的線性模型的參數是有一定的差別的 相對於最小二乘法(OLS)二樣，貝葉斯脊回歸得到的參數比較接近於0 貝葉斯脊回歸的先驗分布是參數向量的高斯分布 """ print(__doc__) import numpy as np import matplotlib.pyplot as plt from scipy import stats import time from sklearn.linear_model import BayesianRidge, LinearRegression ############################################################################### # 隨機函數的種子 np.random.seed(int(time.time()) % 100) # 樣本數目為100，特征數目也是100 n_samples, n_features = 100, 100 # 生成高斯分布 X = np.random.randn(n_samples, n_features) # 首先使用alpha為4的先驗分布. alpha_ = 4. w = np.zeros(n_features) # 隨機提取10個特征出來作為樣本特征 relevant_features = np.random.randint(0, n_features, 10) # 基於先驗分布，產生特征對應的初始權值 for i in relevant_features: w[i] = stats.norm.rvs(loc=0, scale=1. / np.sqrt(alpha_)) # 產生alpha為50的噪聲 alpha_ = 50. noise = stats.norm.rvs(loc=0, scale=1. / np.sqrt(alpha_), size=n_samples) # 產生目標數據 y = np.dot(X, w) + noise ############################################################################### # 使用貝葉斯脊回歸擬合數據 clf = BayesianRidge(compute_score=True) clf.fit(X, y) # 使用最小二乘法擬合數據 ols = LinearRegression() ols.fit(X, y) ############################################################################### # 作圖比較兩個方法的結果 plt.figure(figsize=(6, 5)) plt.title("Weights of the model") plt.plot(clf.coef_, 'b-', label="Bayesian Ridge estimate") plt.plot(w, 'g-', label="Ground truth") plt.plot(ols.coef_, 'r--', label="OLS estimate") plt.xlabel("Features") plt.ylabel("Values of the weights") plt.legend(loc="best", prop=dict(size=12)) plt.show()

其運行結果為：

貝葉斯脊回歸參數：
[ -7.16688614e-02 3.73638195e-02 -4.04171217e-02 7.28338457e-03 ... 2.60774221e-01 4.26079127e-02 1.01660304e-02 6.79853349e-02] 最小二乘法參數： [-1.77270023 -0.38832798 0.58907738 -1.61514115 0.58202424 0.09483505 ... -1.37056305 2.81533169 0.02429617 0.90196961]

可以很容易的看出，相對於最小二乘法(OLS)二樣，貝葉斯脊回歸得到的參數比較接近於0