頻率主義線性回歸概述
線性回歸的頻率主義觀點可能你已經(jīng)學(xué)過了：該模型假定因變量（y）是權(quán)重乘以一組自變量（x）的線性組合。完整的公式還包含一個誤差項以解釋隨機采樣噪聲。如有兩個自變量時，方程為：
模型中，y是因變量，β是權(quán)重（稱為模型參數(shù)），x是自變量的值，ε是表示隨機采樣噪聲的誤差項或變量的影響。
線性回歸是一個簡單的模型，它可以很容易解釋：是截距項，其他權(quán)重β表示增加自變量對因變量的影響。 例如，如果是1.2，那么對于中的每個單位增加，響應(yīng)將增加1.2。
我們可以使用矩陣方程將線性模型推廣到任意數(shù)量的預(yù)測變量。 在預(yù)測矩陣中添加一個常數(shù)項1以解釋截距，我們可以將矩陣公式寫為：
從訓(xùn)練數(shù)據(jù)中學(xué)習(xí)線性模型的目標(biāo)是找到最能解釋數(shù)據(jù)的系數(shù)β。 在頻率主義線性回歸中，最好的解釋是采用殘差平方和（RSS）的系數(shù)β。 RSS是已知值（y）和預(yù)測模型輸出之間的差值的總和（?，表示估計的明顯的y-hat）。 殘差平方和是模型參數(shù)的函數(shù)：
總和被用于訓(xùn)練集中的N個數(shù)據(jù)點。 我們在這里不會詳細討論這個細節(jié)，但是這個方程對于模型參數(shù)β有封閉解，可以使誤差最小化。 這被稱為β的最大似然估計，因為它是給定輸入X和輸出y的最可能的值。 以矩陣形式表示的封閉形式解為：
（再一次，我們必須在β上放上'帽子'，因為它代表了模型參數(shù)的估計值。）不要讓矩陣算術(shù)嚇跑你！ 感謝像Python中的Scikit-learn這樣的庫，我們通常不需要手工計算（盡管編碼線性回歸是一種很好的做法）。 這種通過最小化RSS來擬合模型參數(shù)的方法稱為最小二乘法（OLS）。
我們從頻率主義線性回歸中得到的僅僅是基于訓(xùn)練數(shù)據(jù)的模型參數(shù)的單一估計。 我們的模型完全被數(shù)據(jù)告知：在這個視圖中，我們需要知道的模型的所有信息都編碼在我們可用的訓(xùn)練數(shù)據(jù)中。
一旦我們有了β-hat，我們可以通過應(yīng)用我們的模型方程來估計任何新數(shù)據(jù)點的輸出值：
作為OLS的一個例子，我們可以對真實世界的數(shù)據(jù)進行線性回歸，這些數(shù)據(jù)的持續(xù)時間和消耗的熱量為15000次運動觀察。 以下是通過求解上述模型參數(shù)的矩陣方程得到的數(shù)據(jù)和OLS模型：
使用OLS，我們得到模型參數(shù)的單個估計值，在這種情況下，線的截距和斜率。我們可以寫出由OLS產(chǎn)生的等式：
從斜坡上，我們可以說每一分鐘的鍛煉就能燃燒7.17卡路里。 這種情況下的截距并不有用，因為它告訴我們，如果我們運動0分鐘，我們會燃燒-21.86卡路里！ 這只是OLS擬合程序的一個人為因素，它找到了盡可能減少訓(xùn)練數(shù)據(jù)錯誤的線條，無論它是否物理上合理。
如果我們有一個新的數(shù)據(jù)點，說一個15.5分鐘的運動持續(xù)時間，我們可以將其插入到方程式中，以獲得燃燒卡路里的點估計值：
最小二乘法給出了輸出的單點估計，我們可以將其解釋為給定數(shù)據(jù)的最可能估計。 但是，如果我們有一個小數(shù)據(jù)集，我們可能會將我們的估計值表示為可能值的分布，這就是貝葉斯線性回歸。
完整代碼：
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不被估計為單個值，而是被假定為從正態(tài)分布中抽取。 貝葉斯線性回歸模型是：
輸出y由一個以均值和方差為特征的正態(tài)（高斯）分布產(chǎn)生。 線性回歸的均值是權(quán)重矩陣乘以預(yù)測矩陣的轉(zhuǎn)置。 方差是標(biāo)準(zhǔn)差σ的平方（乘以恒等矩陣，因為這是模型的多維表達式）。
貝葉斯線性回歸的目的不是找到模型參數(shù)的單一“最佳”值，而是確定模型參數(shù)的后驗分布。 不僅是由概率分布產(chǎn)生的響應(yīng)，而且假定模型參數(shù)也來自分布。 模型參數(shù)的后驗概率取決于訓(xùn)練輸入和輸出：
這里，
讓我們停下來思考這意味著什么。與OLS相比，模型參數(shù)有一個后驗分布，它與數(shù)據(jù)乘以參數(shù)先驗概率的可能性成正比。在這里我們可以觀察到貝葉斯線性回歸的兩個主要好處。
1. 先驗：如果我們有領(lǐng)域知識，或者猜測模型參數(shù)應(yīng)該是什么，那么我們可以將它們包括在我們的模型中，這與頻率方法不同，后者假設(shè)所有參數(shù)都來自數(shù)據(jù)。如果我們沒有提前做出任何估計，那么我們可以使用非信息性的先驗來確定正態(tài)分布等參數(shù)。
2. 后驗：執(zhí)行貝葉斯線性回歸的結(jié)果是基于數(shù)據(jù)和先驗的可能模型參數(shù)的分布。這使我們能夠量化我們對模型的不確定性：如果數(shù)據(jù)少，后驗分布將更加分散。
隨著數(shù)據(jù)點數(shù)量的增加，可能性會沖刷先驗，并且在無限數(shù)據(jù)的情況下，參數(shù)的輸出會收斂到從OLS獲得的值。
作為分布的模型參數(shù)的表達形式包含了貝葉斯的世界觀：我們從最初的估計開始，即先驗，并且隨著我們收集更多的證據(jù)，我們的模型變得不那么錯了。貝葉斯推理是我們直覺的自然延伸。通常，我們有一個最初的假設(shè)，當(dāng)我們收集支持或反駁我們想法的數(shù)據(jù)時，我們改變了我們的世界模型（理想情況下這是我們的理由）！
實現(xiàn)貝葉斯線性回歸
在實踐中，評估模型參數(shù)的后驗分布對于連續(xù)變量是難以處理的，所以我們使用抽樣方法從后面抽取樣本以近似后驗。從分布中抽取隨機樣本以近似分布的技術(shù)是蒙特卡羅方法的一種應(yīng)用。有許多用于蒙特卡羅采樣的算法，其中最常見的是馬爾可夫鏈蒙特卡洛變體。
貝葉斯線性建模應(yīng)用
我將跳過本文的代碼，但實現(xiàn)貝葉斯線性回歸的基本過程是：為模型參數(shù)指定先驗（我在本例中使用了正態(tài)分布），創(chuàng)建模型映射訓(xùn)練輸入到訓(xùn)練輸出，然后用馬爾可夫鏈蒙特卡羅（MCMC）算法從后驗分布中抽取樣本作為模型參數(shù)。最終結(jié)果將是參數(shù)的后驗分布。我們可以檢查這些分布以了解發(fā)生了什么。
第一個圖顯示模型參數(shù)的后驗分布的近似值。這些是MCMC 1000步的結(jié)果，這意味著算法從后驗分布中抽取了1000步。
如果我們將斜率和截距的平均值與OLS得到的平均值進行比較（OLS的截距為-21.83，斜率為7.17），會發(fā)現(xiàn)它們非常相似。但是，盡管我們可以將均值用作單點估計，但我們也可以為模型參數(shù)提供一系列可能的值。隨著數(shù)據(jù)點數(shù)量的增加，這個范圍將縮小并收斂一個代表模型參數(shù)更大置信度的值。 （在貝葉斯推斷中，變量的范圍稱為可信區(qū)間，與頻率推理中的置信區(qū)間的解釋略有不同）。
當(dāng)我們想用貝葉斯模型進行線性擬合時，我們可以繪制一系列線條，而不是僅顯示估計值，每條線條表示模型參數(shù)的不同估計值。隨著數(shù)據(jù)點數(shù)量的增加，線條開始重疊，因為模型參數(shù)中的不確定性逐漸減小。
為了證明模型中數(shù)據(jù)點的數(shù)量的影響，我使用了兩個模型，第一個模型，使用了500個數(shù)據(jù)點，第二個使用了15000個數(shù)據(jù)點。每個圖表顯示了100個從參數(shù)分布中抽樣的模型。
當(dāng)使用較少的數(shù)據(jù)點時，擬合中的變化更大，這表示模型中存在更大的不確定性。 有了所有的數(shù)據(jù)點，OLS和貝葉斯擬合幾乎完全相同，因為數(shù)據(jù)的可能性使得先驗數(shù)據(jù)逐漸被覆蓋。
當(dāng)使用我們的貝葉斯線性模型預(yù)測單個數(shù)據(jù)點的輸出時，我們得到的仍是一個分布。 以下是運行15.5分鐘的卡路里數(shù)量的概率密度圖。 紅色垂直線表示來自O(shè)LS的點估計。
結(jié)論
貝葉斯主義和頻率主義的沒有好壞，只有針對它們的合適的應(yīng)用場景。
在數(shù)據(jù)有限或者我們想要在模型中使用某些先驗知識的問題中，貝葉斯線性回歸方法既可以包含先前的信息，也可以表明我們的不確定性。 貝葉斯線性回歸反映了貝葉斯框架：當(dāng)我們收集更多數(shù)據(jù)時，可以改進對數(shù)據(jù)的估計。 貝葉斯觀點是一種直觀的觀察世界的方式，貝葉斯推理可以成為其頻率主義者推論的有用替代方案。 數(shù)據(jù)科學(xué)不是偏袒某一方，而是為了找出工作的最佳工具，并且使用更多的技巧才能使你更有效！
參考鏈接：
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
