2018/11/28

Multiple regression and Polynomial regression (多重回歸與多項式回歸)

05 Multiple regression and Polynomial regression (多重回歸與多項式回歸)

Multiple linear regression (多重線性回歸)

∎ 直接求解
考慮多個features,例如目標的 $y$ 值可寫成
$y=\beta _0+\beta _1x_1+\cdots +\beta _mx_m$,
其中 $\boldsymbol{x}$ 是特徵,$\boldsymbol{\beta}$ 是要求的係數,此式表示共有$m$個特徵。在簡單回歸中我們求$\beta _0$和 $\beta _1$,並用 $y=\beta _0+\beta _1x$ 這直線來預測新成員,其中 $\beta _0$是 intercept (截距)而$\beta _1$ 是直線的斜率,在此情形我們是用直線當系統模型來預測未知的值,且此情況只有一個特徵,直線可畫在平面上,若有兩個特徵則直線就要建立在立體空間裡,更多維的情況依此類推,若有$m$維就必需在$m+1$維的空間畫一條直線,亦即不管有幾個特徵,我們都只用一直線當預測的模型。

以兩個特徵,$n$個實例為例子,表示成
$y_0=\beta _0+\beta _1x_{01}+\beta _2x_{02}$
$y_1=\beta _0+\beta _1x_{11}+\beta _2x_{12}$
$\cdots$
$y_{n-1}=\beta _0+\beta _1x_{n-1,1}+\beta _2x_{n-1,2}$
其中$x_{i,1}$及$x_{i,2}$是第$i$個實例的兩個特徵值,若此例中的$n \leq 3$,則至少可得一解來滿足所有給定的方程式,但在多重回歸中考慮的是$n>3$的情形,亦即我們可能無法找出任何向量 $\boldsymbol{\beta}$ 滿足所有給定的$n$個方程式,只能求出一個向量帶回式子中使得Lost function (一般是平均平方誤差)最小。

若我們把上述寫成矩陣形式,得
$\textbf{Y}=\textbf{X}\boldsymbol\beta$
其中$\textbf{X}$的第一行全為1,維度是$n \times 3$,而$\textbf{Y}$的維度是$n \times 1$,$\boldsymbol\beta$是$3 \times 1$
由線性代數知$\boldsymbol\beta$可由下列求得
$\boldsymbol\beta = (\textbf{X}^T \textbf{X})^{-1}\textbf{X}^T\textbf{Y}$

假設有一組披薩的價格資料賞為訓練集如下,

實例 直徑(吋) 配料 價格
1 6 2 7
2 8 1 9
3 10 0 13
4 14 2 17.5
5 18 0 18

直接套用$\boldsymbol\beta $公式
import numpy as np
from numpy.linalg import inv
from numpy import dot, transpose
X = [[1, 6, 2], [1, 8, 1], [1, 10, 0], [1, 14, 2], [1, 18, 0]]
y = [[7], [9], [13], [17.5], [18]]
inv1 = inv(np.matmul(transpose(X), X))
aa=np.matmul(inv1,transpose(X))
beta = np.matmul(aa,y)
print(beta)
得 [[1.1875    ]
    [1.01041667]
    [0.39583333]]
若我們有一筆資料(7,2),可用已求得的$\boldsymbol\beta $來預測,如
y1=[1,7,2]
predict=np.matmul(y1,beta)
print(predict)
得[9.05208333]
此外, numpy亦有least square function 可以直接求解,
from numpy.linalg import lstsq
print(lstsq(X, y, rcond=None)[0])
得   [[1.1875    ]
       [1.01041667]
       [0.39583333]]
函數lstsq()中的rcond是用於設定浮點數,若不設None會有警告出現,且此函數的輸出有很多貢目,用[0]取出解。

 ∎ 用sklearn求解
若用同上的訓練集以及下列的測試集,

實例 直徑(吋) 配料 價格
1 8 2 11
2 9 0 8.5
3 11 2 15
4 16 2 18
512 0 11


from sklearn.linear_model import LinearRegression
X = [[6, 2], [8, 1], [10, 0], [14, 2], [18, 0]]
y = [[7], [9], [13], [17.5], [18]]
model = LinearRegression()
model.fit(X, y)
print(model.intercept_)
print(model.coef_)
X_test = [[8, 2], [9, 0], [11, 2], [16, 2], [12, 0]]
y_test = [[11], [8.5], [15], [18], [11]]
predictions = model.predict(X_test)
for i, prediction in enumerate(predictions):
   print('Predicted: %s, Target: %s' % (prediction, y_test[i]))
print('R-squared: %.2f' % model.score(X_test, y_test))
得 [1.1875]
     [[1.01041667 0.39583333]]
     Predicted: [10.0625], Target: [11]
     Predicted: [10.28125], Target: [8.5]
     Predicted: [13.09375], Target: [15]
     Predicted: [18.14583333], Target: [18]
     Predicted: [13.3125], Target: [11]
     R-squared: 0.77
可見R-squared 由之前 的 0.6620提升到0.77,表示增加feature對預測的結果有改善。其中使用了LinearRegression()為model,用fit去訓練 ,用predict預測,最後用score列印R-squared。

Polynomial regression (多項式回歸)

∎ 一個特徵二階多項式
線性回歸對每一個特徵使用一個變數,利用直線來預測,而在多項式回歸則是用曲線來預測 ,曲線的複雜度和特徵的數目使用多項式的階次有關。例如一個特徵二階多項式就是建立下列曲線用以預測,
$y=\beta_0+\beta_1x+\beta_2x^2$
以下列的資料集來比較簡單回歸和多項式回歸,
訓練集如下,
實例 直徑(吋) 價格
1 6 7
2 8 9
3 10 13
4 14 17.5
5 18 18

測試集如下

實例 直徑(吋) 價格
1 6 8
2 8 12
3 11 15
4 16 18

下列是用以比較的程式
import numpy as np
import matplotlib.pyplot as plt
from sklearn.linear_model import LinearRegression
from sklearn.preprocessing import PolynomialFeatures
X_train = [[6], [8], [10], [14], [18]]
y_train = [[7], [9], [13], [17.5], [18]]
X_test = [[6], [8], [11], [16]]
y_test = [[8], [12], [15], [18]]
regressor = LinearRegression()
regressor.fit(X_train, y_train)
xx = np.linspace(0, 26, 100)
yy = regressor.predict(xx.reshape(xx.shape[0], 1))
plt.plot(xx, yy)
quadratic_featurizer = PolynomialFeatures(degree=2)
X_train_quadratic = quadratic_featurizer.fit_transform(X_train)
X_test_quadratic = quadratic_featurizer.transform(X_test)
regressor_quadratic = LinearRegression()
regressor_quadratic.fit(X_train_quadratic, y_train)
xx_quadratic = quadratic_featurizer.transform(xx.reshape(xx.shape[0], 1))
plt.plot(xx, regressor_quadratic.predict(xx_quadratic), c='r', linestyle='--')
plt.title('Pizza price regressed on diameter')
plt.xlabel('Diameter in inches')
plt.ylabel('Price in dollars')
plt.axis([0, 25, 0, 25])
plt.grid(True)
plt.scatter(X_train, y_train)
plt.show()
print('Simple linear regression R-squared', regressor.score(X_test, y_test))
print('Quadratic regression R-squared',
     regressor_quadratic.score(X_test_quadratic, y_test))

Simple linear regression R-squared 0.809726797707665
Quadratic regression R-squared 0.8675443656345073

可見多項式回歸的R-squared比簡單線性回歸好。
∎ PolynomialFeatures()函數
程式中的 PolynomialFeatures()是用來把輸入的data建構成多項式的形式,再用於多特徵的線性回歸,以完成多項式回歸的工作。例如若有$a$,$b$兩個特徵,那麼2次多項式為$1,a,b, a^2 ,ab,  b^2$。 此例中只有一個特徵,是先把6轉成[[ 1., 6., 36.], 把8 轉成[ 1., 8., 64.],.....。此函數有三個參數:
  • degree:控制多項式的階次(order)
  • interaction_only:預設為False,如果指定為True,那麼就不會有特徵自己和自己結合的項,上面例子的二次項中沒有 $a^2$ 和 $b^2$。
  • include_bias:預設為True。如果為False的話,上面例子中就不會有1那一項。

∎ Overfitting (過適配)
上例若把多項式的階次改成9,執行後得
degree=9
Simple linear regression R-squared 0.809726797707665
Quadratic regression R-squared -0.09435666704293633

model 太複雜,系統記住training data但對預測資料的性能太差,上圖中R-squared 只有-0.09,對16吋預測值低於10很顯然不合理

∎ Regularization(調適)
用於解決 overfitting 的問題,對複雜度加入懲罰項目使得太複雜的模型的Cost function反而變差,以常用的Ridge regression為例,Cost function定義如下:
$RSS_{\rm{ridge}}=\sum_{i=0}^{n-1}(y_i-x_i\beta_i)^2+\lambda \sum_{j=0}^{p}\beta_i^2$
sklearn裡叫用的例子
from sklearn.linear_model import Ridge  
from sklearn.linear_model import LinearRegression
import numpy as np  
n_samples, n_features = 10, 1  
np.random.seed(0)  
y = np.random.randn(n_samples)  
X = np.random.randn(n_samples, n_features)  
regressor = LinearRegression()
regressor.fit(X, y)
clf = Ridge(alpha=1)  
clf.fit(X, y)
print(clf.intercept_)
print(clf.coef_)
print(regressor.intercept_)
print(regressor.coef_)
0.7042107630757627
[0.08439472]
0.6967543540472676
[0.10300568]

Gradient descent (梯度下降) 演算法

梯度下降法是一個一階最佳化算法,用以尋找函數的極小值。要尋找函數的區域極小值(local minimum),必須向函數上當前點對應梯度(或者是近似梯度)的反方向(正方向為上升方向),以規定步長躍進,到新點之後再以相同的方法不斷躍進,經不斷迭代後若梯度小於一事先定義的數值則迭代停止,表示已約略得到區域極小值。因Cost function是convex function (凸函數)所以區域極小值也是整體的極小值。

例如函數$F(\textbf x)$是個多維變數函數,假設在點$a$鄰近可微分,在此函數於$a$點向負梯度的方向遞減最快速,所以若以迭代的方式,令
$\textbf {a}_{n+1}=\textbf {a}_n-\gamma_n\Delta F(\textbf {a}_n)$

則$F(\textbf {a}_n) \geq  F(\textbf {a}_{n+1})$,其中$\gamma_n$為一小的正實數,可以是常數或隨$n$而變。當以一任意選擇的啟始點開始迭代,歷經多次迭代以後若遞減的量小於一定數值就停止迭,亦即已找到理想的最小值。在機器學習的過程中,要最佳化的函數 $F(\textbf x)$ 中的 $\textbf x$ 並不是連續空間,而是限定在所給的樣本(實例)的離散點,Gradient descent 在機器學習中的應用以如何取用所給的樣本求取梯度用以進行迭代,可分成三大類:

  • Batch gradient descent (BGD, 整批梯度下降):每次迭代時所有的樣本都參與更新的運算,若樣本數目多時運算速度會很慢。
  • Stochastic gradient descent (SGD, 隨機梯度下降):每次迭代時隨機選取一個樣本進行計算及參數更新 ,若樣本數目多時此法較適合。
  • Mini-batch gradient descent (MBGD, 小批量梯度下降):每次迭代取一定的數量(如10個或20個,但不是全部樣本)參與更新運算,此法為上兩法的折衷。

梯度下降演算法在機器學習及深度學習扮演很重要角色,後續會用得到。


參考文獻
Gavin Hackeling, Mastering Machine Learning with scikit-learn, 2nd, Packt Publishing, 2017
https://en.wikipedia.org/wiki/Gradient_descent

沒有留言:

張貼留言