貝葉斯線性迴歸簡介（附完整程式碼）

今晚（4月25日）8點，七月線上公開課【如何從零轉崗AI】，點選文末“閱讀原文”進入直播間。

作者 | William Koehrsen

編譯 | 專知

參與 | Yingying, Xiaowen

【導讀】應用貝葉斯推理的重點領域之一是貝葉斯線性模型。我們首先簡要回顧一下頻率主義學派的線性迴歸方法，接著介紹貝葉斯推斷，並試著應用於簡單的資料集。

Introduction to Bayesian Linear Regression

頻率主義線性迴歸概述

---

---

 

線性迴歸的頻率主義觀點可能你已經學過了：該模型假定因變數（y）是權重乘以一組自變數（x）的線性組合。完整的公式還包含一個誤差項以解釋隨機取樣噪聲。如有兩個自變數時，方程為：

模型中，y是因變數，β是權重（稱為模型引數），x是自變數的值，ε是表示隨機取樣噪聲的誤差項或變數的影響。

線性迴歸是一個簡單的模型，它可以很容易解釋：是截距項，其他權重β表示增加自變數對因變數的影響。 例如，如果是1.2，那麼對於中的每個單位增加，響應將增加1.2。

我們可以使用矩陣方程將線性模型推廣到任意數量的預測變數。 在預測矩陣中新增一個常數項1以解釋截距，我們可以將矩陣公式寫為：

從訓練資料中學習線性模型的目標是找到最能解釋資料的係數β。 在頻率主義線性迴歸中，最好的解釋是採用殘差平方和（RSS）的係數β。 RSS是已知值（y）和預測模型輸出之間的差值的總和（ŷ，表示估計的明顯的y-hat）。 殘差平方和是模型引數的函式：

總和被用於訓練集中的N個資料點。 我們在這裡不會詳細討論這個細節，但是這個方程對於模型引數β有封閉解，可以使誤差最小化。 這被稱為β的最大似然估計，因為它是給定輸入X和輸出y的最可能的值。 以矩陣形式表示的封閉形式解為：

（再一次，我們必須在β上放上'帽子'，因為它代表了模型引數的估計值。）不要讓矩陣算術嚇跑你！ 感謝像Python中的Scikit-learn這樣的庫，我們通常不需要手工計算（儘管編碼線性迴歸是一種很好的做法）。 這種通過最小化RSS來擬合模型引數的方法稱為最小二乘法（OLS）。

我們從頻率主義線性迴歸中得到的僅僅是基於訓練資料的模型引數的單一估計。 我們的模型完全被資料告知：在這個檢視中，我們需要知道的模型的所有資訊都編碼在我們可用的訓練資料中。

一旦我們有了β-hat，我們可以通過應用我們的模型方程來估計任何新資料點的輸出值：

作為OLS的一個例子，我們可以對真實世界的資料進行線性迴歸，這些資料的持續時間和消耗的熱量為15000次運動觀察。 以下是通過求解上述模型引數的矩陣方程得到的資料和OLS模型：

使用OLS，我們得到模型引數的單個估計值，在這種情況下，線的截距和斜率。我們可以寫出由OLS產生的等式：

從斜坡上，我們可以說每一分鐘的鍛鍊就能燃燒7.17卡路里。 這種情況下的截距並不有用，因為它告訴我們，如果我們運動0分鐘，我們會燃燒-21.86卡路里！ 這只是OLS擬合程式的一個人為因素，它找到了儘可能減少訓練資料錯誤的線條，無論它是否物理上合理。

如果我們有一個新的資料點，說一個15.5分鐘的運動持續時間，我們可以將其插入到方程式中，以獲得燃燒卡路里的點估計值：

最小二乘法給出了輸出的單點估計，我們可以將其解釋為給定資料的最可能估計。 但是，如果我們有一個小資料集，我們可能會將我們的估計值表示為可能值的分佈，這就是貝葉斯線性迴歸。

完整程式碼：

1. Load in Exercise Data

exercise = pd.read_csv('data/exercise.csv')
calories = pd.read_csv('data/calories.csv')
df = pd.merge(exercise, calories, on = 'User_ID')
df = df[df['Calories'] < 300]
df = df.reset_index()
df['Intercept'] = 1
df.head()

2. Plot Relationship

plt.figure(figsize=(8, 8))

plt.plot(df['Duration'], df['Calories'], 'bo');
plt.xlabel('Duration (min)', size = 18); plt.ylabel('Calories', 
size = 18); 
plt.title('Calories burned vs Duration of Exercise', size = 20);

# Create the features and response
X = df.loc[:, ['Intercept', 'Duration']]
y = df.ix[:, 'Calories']

3. Implement Ordinary Least Squares Linear Regression by Hand

# Takes a matrix of features (with intercept as first column) 
# and response vector and calculates linear regression coefficients
def linear_regression(X, y):
# Equation for linear regression coefficients
    beta = np.matmul(np.matmul(np.linalg.inv(np.matmul(X.T, X)), X.T),
y)
return beta

# Run the by hand implementation
by_hand_coefs = linear_regression(X, y)
print('Intercept calculated by hand:', by_hand_coefs[0])
print('Slope calculated by hand: ', by_hand_coefs[1])

xs = np.linspace(4, 31, 1000)
ys = by_hand_coefs[0] + by_hand_coefs[1] * xs

plt.figure(figsize=(8, 8))
plt.plot(df['Duration'], df['Calories'], 'bo', label = 'observations', 
alpha = 0.8);
plt.xlabel('Duration (min)', size = 18); plt.ylabel('Calories', size = 18); 
plt.plot(xs, ys, 'r--', label = 'OLS Fit', linewidth = 3)
plt.legend(prop={'size': 16})
plt.title('Calories burned vs Duration of Exercise', size = 20);

貝葉斯線性迴歸

---

---

 

在貝葉斯觀點中，我們使用概率分佈而不是點估計來進行線性迴歸。 y不被估計為單個值，而是被假定為從正態分佈中抽取。 貝葉斯線性迴歸模型是：

輸出y由一個以均值和方差為特徵的正態（高斯）分佈產生。 線性迴歸的均值是權重矩陣乘以預測矩陣的轉置。 方差是標準差σ的平方（乘以恆等矩陣，因為這是模型的多維表示式）。

貝葉斯線性迴歸的目的不是找到模型引數的單一“最佳”值，而是確定模型引數的後驗分佈。 不僅是由概率分佈產生的響應，而且假定模型引數也來自分佈。 模型引數的後驗概率取決於訓練輸入和輸出：

這裡，是給定輸入和輸出的模型引數的後驗概率分佈。 這等於輸出乘以給定輸入的β的先驗概率併除以歸一化常數的可能性。 這是貝葉斯定理的一個簡單表示式，貝葉斯定理是貝葉斯推理的基本基礎：

讓我們停下來思考這意味著什麼。與OLS相比，模型引數有一個後驗分佈，它與資料乘以引數先驗概率的可能性成正比。在這裡我們可以觀察到貝葉斯線性迴歸的兩個主要好處。

1. 先驗：如果我們有領域知識，或者猜測模型引數應該是什麼，那麼我們可以將它們包括在我們的模型中，這與頻率方法不同，後者假設所有引數都來自資料。如果我們沒有提前做出任何估計，那麼我們可以使用非資訊性的先驗來確定正態分佈等引數。

2. 後驗：執行貝葉斯線性迴歸的結果是基於資料和先驗的可能模型引數的分佈。這使我們能夠量化我們對模型的不確定性：如果資料少，後驗分佈將更加分散。

隨著資料點數量的增加，可能性會沖刷先驗，並且在無限資料的情況下，引數的輸出會收斂到從OLS獲得的值。

作為分佈的模型引數的表達形式包含了貝葉斯的世界觀：我們從最初的估計開始，即先驗，並且隨著我們收集更多的證據，我們的模型變得不那麼錯了。貝葉斯推理是我們直覺的自然延伸。通常，我們有一個最初的假設，當我們收集支援或反駁我們想法的資料時，我們改變了我們的世界模型（理想情況下這是我們的理由）！

實現貝葉斯線性迴歸

---

---

 

在實踐中，評估模型引數的後驗分佈對於連續變數是難以處理的，所以我們使用抽樣方法從後面抽取樣本以近似後驗。從分佈中抽取隨機樣本以近似分佈的技術是蒙特卡羅方法的一種應用。有許多用於蒙特卡羅取樣的演算法，其中最常見的是馬爾可夫鏈蒙特卡洛變體。

貝葉斯線性建模應用

---

---

 

我將跳過本文的程式碼，但實現貝葉斯線性迴歸的基本過程是：為模型引數指定先驗（我在本例中使用了正態分佈），建立模型對映訓練輸入到訓練輸出，然後用馬爾可夫鏈蒙特卡羅（MCMC）演算法從後驗分佈中抽取樣本作為模型引數。最終結果將是引數的後驗分佈。我們可以檢查這些分佈以瞭解發生了什麼。

第一個圖顯示模型引數的後驗分佈的近似值。這些是MCMC 1000步的結果，這意味著演算法從後驗分佈中抽取了1000步。

如果我們將斜率和截距的平均值與OLS得到的平均值進行比較（OLS的截距為-21.83，斜率為7.17），會發現它們非常相似。但是，儘管我們可以將均值用作單點估計，但我們也可以為模型引數提供一系列可能的值。隨著資料點數量的增加，這個範圍將縮小並收斂一個代表模型引數更大置信度的值。 （在貝葉斯推斷中，變數的範圍稱為可信區間，與頻率推理中的置信區間的解釋略有不同）。

當我們想用貝葉斯模型進行線性擬合時，我們可以繪製一系列線條，而不是僅顯示估計值，每條線條表示模型引數的不同估計值。隨著資料點數量的增加，線條開始重疊，因為模型引數中的不確定性逐漸減小。

為了證明模型中資料點的數量的影響，我使用了兩個模型，第一個模型，使用了500個資料點，第二個使用了15000個資料點。每個圖表顯示了100個從引數分佈中抽樣的模型。 

當使用較少的資料點時，擬閤中的變化更大，這表示模型中存在更大的不確定性。 有了所有的資料點，OLS和貝葉斯擬合幾乎完全相同，因為資料的可能性使得先驗資料逐漸被覆蓋。

當使用我們的貝葉斯線性模型預測單個資料點的輸出時，我們得到的仍是一個分佈。 以下是執行15.5分鐘的卡路里數量的概率密度圖。 紅色垂直線表示來自OLS的點估計。

 

結論

---

---

 

貝葉斯主義和頻率主義的沒有好壞，只有針對它們的合適的應用場景。

在資料有限或者我們想要在模型中使用某些先驗知識的問題中，貝葉斯線性迴歸方法既可以包含先前的資訊，也可以表明我們的不確定性。 貝葉斯線性迴歸反映了貝葉斯框架：當我們收集更多資料時，可以改進對資料的估計。 貝葉斯觀點是一種直觀的觀察世界的方式，貝葉斯推理可以成為其頻率主義者推論的有用替代方案。 資料科學不是偏袒某一方，而是為了找出工作的最佳工具，並且使用更多的技巧才能使你更有效！

參考連結：

1.https://www.quantstart.com/articles/Bayesian-Linear-Regression-Models-with-PyMC3

2.http://twiecki.github.io/blog/2013/08/12/bayesian-glms-1/

3.https://wiseodd.github.io/techblog/2017/01/05/bayesian-regression/

4.PyMC3 Introduction

完整程式碼連結：

https://github.com/WillKoehrsen/Data-Analysis/blob/master/bayesian_lr/Bayesian%20Linear%20Regression%20Demonstration.ipynb

原文連結：

https://towardsdatascience.com/introduction-to-bayesian-linear-regression-e66e60791ea7

★免費公開課★

今晚8點，掃碼進入直播間