WayneLinear regression 是一個資料分析技術，且使用線性函數來預測未知的資料。雖然 linear regression 模型相對簡單，但它是一個成熟的統計技術。
Table of Contents
線性回歸
線性回歸模型（Linear Regression Model）
Linear regression 的 model function 如下。其中 w 和 b 是參數。當我們將變數 x 帶入後，函數 fw,b 會回傳一個預測值 ŷ（y hat）。
要注意的是，fw,b 回傳的是一個預測值 ŷ（prediction value），而非真正的值 y（target value），如下圖所示。當 fw,b 預測出來的 ŷ 很接近 y 時，則我們可以說 fw,b 的準確率很高。我們透過調整 w 和 b 來提高 fw,b 的準確率。
以下的程式碼中，f_wb() 實作了 linear regression model。
def f_wb(x, w, b):
"""
Linear regression model
Parameters
----------
x: (ndarray (m,)) or (scalar)
m is the number of samples
w: (scalar)
b: (scalar)
Returns
-------
(ndarray (m,)) or (scalar)
"""
return w * x + b成本函數（Cost Function）
在介紹如何求得 w 和 b 之前，我們要先了解什麼是 cost function。Cost function 是用來衡量 fw,b 的 prediction values 與 target values 之間的殘差（residual），也就是它們之間的差值。有了 cost function，我們就可以衡量調整後的 w 和 b 是否比原來的好。
Linear regression 的 cost function 是 squared error，如下。
我們最終想要找到一組 w 和 b，使得 J(w,b) 會是最小。
以下的程式碼中，compute_cost() 實作了 cost function J(w,b)。
import numpy as np
def compute_cost(x, y, w, b):
"""
Compute cost
Parameters
----------
x: (ndarray (m,))
m is the number of samples
y: (ndarray (m,))
w: (scalar)
b: (scalar)
Returns
-------
(scalar)
"""
m = x.shape[0]
y_hat = f_wb(x, w, b)
return 1 / (2 * m) * np.sum((y_hat - y) ** 2)求解參數 w 和 b
Linear regression 的 model function 很簡單，然而問題是要如何求解 w 和 b，並使得 J(w,b) 會是最小呢？在 machine learning 中，我們使用的是梯度下降（gradient descent）來求解 w 和 b。不過我們會先介紹另外兩種方法，來幫助讀者了解 gradient descent 的優點。
在介紹各個方法時，我們會使用以下的範例來講解。
import matplotlib.pyplot as plt
import numpy as np
x_train = np.array([10, 20, 32, 44, 66])
y_train = np.array([10, 30, 40, 44, 54])
plt.scatter(x_train, y_train, 60, marker='o', c='c', label='y_train')
plt.xlabel('x_train')
plt.ylabel('y_train')
plt.legend()
聯立方程式（Simultaneous Equations）
我們都知道平面上的任意兩點可構成一條直線。我們挑選 (10, 10) 和 (66, 54) 來計算 w 和 b。
我們將計算出來的 w 和 b 代入 model function，並在座標圖上畫出用 simultaneous equations 算出的 model function 以及 prediction values。
# simultaneous equations
w_se = 11 / 14 # 0.7857
b_se = 30 / 14 # 2.1428
y_hat_se = f_wb(x_train, w_se, b_se)
plt.scatter(x_train, y_train, 60, marker='o', c='c', label='y_train')
plt.plot(x_train, y_hat_se, 'ro-', label='y_hat_se')
plt.xlabel('x_train')
plt.ylabel('y_train and y_hat_se')
plt.legend()
以下顯示該 model 的 cost。利用範例中的 5 個點，我們可以計算出 10 組 w 和 b，然後計算它們的 cost。最後，選擇 cost 最小的 w 和 b。
cost_se = compute_cost(x_train, y_train, w_se, b_se) cost_se # Output: # np.float64(36.21836734693877)
使用聯立方程式來求解 w 和 b 看似簡單，然而當 training examples 多的時候，計算會趨於複雜。此外，如果最佳的 w 和 b 的線不會經過任何一個 training example，那我們將無法找到此 w 和 b。
最小平方法（Least Squares）
Least squares 是一個統計學的優化建模方法。它求解的 w 和 b 可使 residual 的平方和最小，也就是 cost 最小。Least squares 的概念是在 model function 中加入一個誤差值 ε（epsilon），來使得 prediction value 可以等於 target value。而且，盡量地使 ε 最小。
將所有的 training examples 代入 model function。可以看出，每一個 training example 的 ε 都可能不同。
接著，計算每個 ε 的平方。
計算所有 ε 的平方的和。然後，我們希望求解出一組 w 和 b，使得所有 ε 的平方的和為最小，也就是越接近零越好。
利用配方法（completing the square）來化簡式子。
我們希望誤差可以最小的話，w 的值應該要如下：
這樣子的話，前面的部分就會等於零。我們繼續用配方法化簡剩下的部分。
要讓上面的式子最小的話，那 b 的值應該要如下，這樣前面的部分就會等於零。
將 b 值代入到 w 的式子裡，並得出 w 的值。
將算出來的 w 和 b 代入 model function，則可以得出我們要的 model。
以下的程式碼中，我們在座標圖上畫出用 least squares 算出的 model function 以及 prediction values。
# least squares
w_ls = 0.788
b_ls = 7.674
y_hat_ls = f_wb(x_train, w_ls, b_ls)
plt.scatter(x_train, y_train, 60, marker='o', c='c', label='y_train')
plt.plot(x_train, y_hat_se, 'm^-', label='y_hat_se')
plt.plot(x_train, y_hat_ls, 'gs-', label='y_hat_ls')
plt.xlabel('x_train')
plt.ylabel('y_train, y_hat_se and y_hat_ls')
plt.legend()
以下顯示該 model 的 cost。
cost_ls = compute_cost(x_train, y_train, w_ls, b_ls) cost_ls # Output # np.float64(15.953221200000002)
相較於用聯立方程式來計算 w 和 b，最小平方法可以算出較好的 w 和 b。但是，當 training examples 很多的時候，整個計算過程會異常複雜。
梯度下降（Gradient Descent）
Gradient descent 是一個最佳化演算法，用來找到一個函數的局部值。其基本原理是，首先，隨便挑選一組 w 和 b，所以可能會在 cost function 的左邊或是右邊，然後下一組選擇谷底方向的 w 和 b。重複此動作直到 J(w,b) 無法再更小。
下一個問題是要如何知道 cost function 的谷底是在哪個方向呢？我們可以對 cost function J(w,b) 做微分，然後就可以得到目前的斜率（slope）。有了斜率，我們就知道要往哪邊移動。
以下為 gradient descent 的演算法。首先，先隨機選一組 w 和 b，或直接設為零。Cost function 的導數乘上一個 learning rate 
會是 w 和 b 要下降的梯度。此外，我們還會給定一個 iteration 次數，代表執行 gradient descent 的次數。因此，我們也許不會得到最佳的 w 和 b。
J(w,b) 對 w 的偏導數（partial derivative）為：
J(w,b) 對 b 的偏導數為：
以下的程式碼實作 J(w,b) 對 w 和 b 的偏導數計算。
def compute_gradient(x, y, w, b):
"""
Compute the gradient for linear regression
Parameters
----------
x: (ndarray (m,))
m is the number of samples
y: (ndarray (m,))
w: (scalar)
b: (scalar)
Returns
-------
dj_dw: (scalar)
dj_db: (scalar)
"""
m = x.shape[0]
dj_dw = 1 / m * np.sum((f_wb(x, w, b) - y) * x)
dj_db = 1 / m * np.sum(f_wb(x, w, b) - y)
return dj_dw, dj_db以下的程式碼實作 gradient descent。其中參數 alpha 是 learning rate，而 epochs 是 iteration 次數。
def perform_gradient_descent(x, y, w_init, b_init, alpha, epochs):
"""
Perform gradient descent
Parameters
----------
x: (ndarray (m,))
m is the number of samples
y: (ndarray (m,))
w_init: (scalar)
b_init: (scalar)
alpha: (float)
epochs: (int)
Returns
-------
w: (scalar)
b: (scalar)
"""
w = w_init
b = b_init
for i in range(epochs):
dj_dw, dj_db = compute_gradient(x, y, w, b)
w -= alpha * dj_dw
b -= alpha * dj_db
return w, b以下的程式碼中，我們令 learning rate 為 0.001，以及 iteration 次數為 3000 次。然後，在座標圖上畫出用 gradient descent 算出的 model function 以及 prediction values。
w_dg, b_dg = perform_gradient_descent(x_train, y_train, 0, 0, 0.001, 3000)
y_hat_gd = f_wb(x_train, w_dg, b_dg)
plt.scatter(x_train, y_train, 60, marker='o', c='c', label='y_train')
plt.plot(x_train, y_hat_se, 'm^-', label='y_hat_se')
plt.plot(x_train, y_hat_ls, 'gs-', label='y_hat_ls')
plt.plot(x_train, y_hat_gd, 'y*-', label='y_hat_gd')
plt.xlabel('x_train')
plt.ylabel('y_train, y_hat_se, y_hat_ls, and y_hat_gd')
plt.legend()
以下顯示該 model 的 cost。
cost_gd = compute_cost(x_train, y_train, w_dg, b_dg) cost_gd # Output # np.float64(18.014426543403744)
雖然在推導 cost function 的導數過程有些複雜，但是整體演算法不會太複雜。當 training examples 多的時候，程式碼也不會變複雜。
以下列出由以上三個方法算出的一些數值，供讀者比較和參考。
(w_se, b_se), (w_ls, b_ls), (w_dg, b_dg) # Output ((0.7857142857142857, 2.142857142857143), (0.788, 7.674), (np.float64(0.8312661524754889), np.float64(5.714916639215474))) (y_train, y_hat_se, y_hat_ls, y_hat_gd) # Output (array([10, 30, 40, 44, 54]), array([10. , 17.85714286, 27.28571429, 36.71428571, 54. ]), array([15.554, 23.434, 32.89 , 42.346, 59.682]), array([14.02757816, 22.34023969, 32.31543352, 42.29062735, 60.5784827 ])) (cost_se, cost_ls, cost_gd) # Output (np.float64(36.21836734693877), np.float64(15.953221200000002), np.float64(18.014426543403744))
Learning Rate
在 gradient descent 中，有一個 learning rate 叫 。由 gradient descent 演算法可以看出，當 較小時，需要較多的次數才能接近谷底，因此效能會較差。反之，則會比較快速地到達底谷。但是，當 太大時，則有可能會無法到達谷底。
多元線性迴歸（Multiple Linear Regression）
至目前為止，我們介紹的是 simple linear regression，只有一個變數，實用性可能不大。接下來，我們將介紹 multiple linear regression。
多元線性迴歸模型（Multiple Linear Regression Model）
Multiple linear regression 有一個以上的變數，如下是 4 個變數的 model function。
所以，multiple linear regression 的 model function 如下。
向量化的 model function 則如下，其中 
是 vectors。
![f_{\vec{w},b}(x^{(i)})=\displaystyle\sum_{j=1}^nw_jx_j+b=\vec{w}\cdot \vec{x}^{(i)}+b \\\\ \vec{w}=[w_1,w_2,\cdots,w_n] \\\\ \vec{x^{(i)}}=[x_1^{(i)},x_2^{(i)},\cdots,x_n^{(i)}]](https://s0.wp.com/latex.php?latex=f_%7B%5Cvec%7Bw%7D%2Cb%7D%28x%5E%7B%28i%29%7D%29%3D%5Cdisplaystyle%5Csum_%7Bj%3D1%7D%5Enw_jx_j%2Bb%3D%5Cvec%7Bw%7D%5Ccdot+%5Cvec%7Bx%7D%5E%7B%28i%29%7D%2Bb+%5C%5C%5C%5C+%5Cvec%7Bw%7D%3D%5Bw_1%2Cw_2%2C%5Ccdots%2Cw_n%5D+%5C%5C%5C%5C+%5Cvec%7Bx%5E%7B%28i%29%7D%7D%3D%5Bx_1%5E%7B%28i%29%7D%2Cx_2%5E%7B%28i%29%7D%2C%5Ccdots%2Cx_n%5E%7B%28i%29%7D%5D+&bg=ffffff&fg=000&s=1&c=20201002)
如果我們將所有的 training examples 放在一起，就會形成一個陣列（matrix）。一般我們會用大寫的 X 來代表 training examples 形成的 matrix。因此，model function 可寫成如下：
向量化的 fw,b 和 simple linear regression 的 model function 很像，它們的實作也差不多。差別在於，向量化的 fw,b 使用內積（dot product）來相乘 X 和 w。
import numpy as np
def f_wb(X, w, b):
"""
Linear regression model
Parameters
----------
X: (ndarray (m, n))
m is the number of samples, n is the number of features
w: (ndarray (n,))
b: (scalar)
Returns
-------
(ndarray (m,))
"""
return np.dot(X, w) + b成本函數（Cost Function）
Multiple linear regression 的 cost function 如下。
以下程式碼中，compute_cost() 實作了向量化的 cost function。
import numpy as np
def compute_cost(X, y, w, b):
"""
Compute cost
Parameters
----------
X: (ndarray (m, n))
m is the number of samples, n is the number of features
y: (ndarray (m,))
w: (ndarray (n,))
b: (scalar)
Returns
-------
(scalar)
"""
m = X.shape[0]
y_hat = f_wb(X, w, b)
return 1 / (2 * m) * np.sum((y_hat - y) ** 2)梯度下降（Gradient Descent）
Multiple linear regression 的 gradient descent 如下。從式子中可以了解到，我們要計算每一個 feature 的參數 wj。然後，計算 J(w,b) 對每一個 wj 的偏導數。
J(w,b) 對每一個 wj 和 b 的偏導數如下：
以下的程式碼實作 J(w,b) 對 wj 和 b 的偏導數計算。
import numpy as np
def compute_gradient(X, y, w, b):
"""
Compute the gradient for linear regression
Parameters
----------
X: (ndarray (m, n))
m is the number of samples, n is the number of features
y: (ndarray (m,))
w: (ndarray (n,))
b: (scalar)
Returns
-------
dj_dw: (ndarray (n,))
dj_db: (scalar)
"""
m = X.shape[0]
y_hat = f_wb(X, w, b)
dj_dw = 1 / m * np.dot(X.T, y_hat - y)
dj_db = 1 / m * np.sum(y_hat - y)
return dj_dw, dj_db以下的程式碼實作 gradient descent。
import numpy as np
def perform_gradient_descent(X, y, w_init, b_init, alpha, epochs):
"""
Perform gradient descent
Parameters
----------
X: (ndarray (m, n))
m is the number of samples, n is the number of features
y: (ndarray (m,))
w_init: (ndarray (n,))
b_init: (scalar)
alpha: (scalar)
epochs: (int)
Returns
-------
w: (ndarray (n,))
b: (scalar)
"""
w = w_init
b = b_init
for i in range(epochs):
dj_dw, dj_db = compute_gradient(X, y, w, b)
w -= alpha * dj_dw
b -= alpha * dj_db
print(f'Epoch {i + 1}, Cost: {compute_cost(X, y, w, b)}')
return w, b以下的程式碼中，我們令 learning rate 為 0.0000001，以及 iteration 次數為 1000 次。然後，輸出 cost 在每個 iteration 的變化。可以看出，cost 在一開始時下降很快，之後越來越慢。
X_train = np.array([[2104, 5, 1, 45], [1416, 3, 2, 40], [852, 2, 1, 35]]) y_train = np.array([460, 232, 178]) w, b = perform_gradient_descent(X_train, y_train, np.zeros(X_train.shape[1]), 0., 0.0000001, 1000) (w, b) # Output # Epoch 1, Cost: 28989.111501449268 # Epoch 2, Cost: 17092.48916583246 # Epoch 3, Cost: 10198.319962562276 # Epoch 4, Cost: 6203.104159134617 # Epoch 5, Cost: 3887.8502848154976 # Epoch 6, Cost: 2546.145030983679 # Epoch 7, Cost: 1768.6173933475704 # Epoch 8, Cost: 1318.0342718091906 # Epoch 9, Cost: 1056.9175763423789 # Epoch 10, Cost: 905.59787958039 # Epoch 11, Cost: 817.9062417617257 # Epoch 12, Cost: 767.0874687152398 # Epoch 13, Cost: 737.6367552876493 # Epoch 14, Cost: 720.5689683005712 # Epoch 15, Cost: 710.6771659424062 # ... # Epoch 995, Cost: 694.9420356037724 # Epoch 996, Cost: 694.9399054268631 # Epoch 997, Cost: 694.9377752873567 # Epoch 998, Cost: 694.9356451852516 # Epoch 999, Cost: 694.9335151205503 # Epoch 1000, Cost: 694.9313850932496 # (array([ 0.20253263, 0.00112386, -0.00213202, -0.00933401]), np.float64(-0.0003572563919114557))
以下顯示利用算出來的 w 和 b 來預測 ŷ，並計算 ŷ 和 y 的 residual。
y_hat = f_wb(X_train, w, b) (y_train, y_hat) # Output # (array([460, 232, 178]), array([425.71175692, 286.41159657, 172.23087045])) cost = compute_cost(X_train, y_train, w, b) cost # Output # np.float64(694.9313850932496)
結語
Linear regression 是一個相較簡單的模型。很適合用來講解模型背後的數學原理。不過，它還是很實用的。Simple linear regression 只有一個變數，實用性可能不大。實際上，我們會使用向量化的 multiple linear regression，因為通常變數會大於一個，而且可能會很多。
參考
- Andrew Ng, Machine Learning Specialization, Coursera.
- 西内啓，機器學習的數學基礎 : AI、深度學習打底必讀，旗標。
Wayne
