多項式回歸
目錄
這篇學習筆記記錄一下由線性模型擴展至非線性模型的多項式回歸。
線性回歸模型形式簡單，有很好的解釋性，但它有不少假設前提，其中最重要的一條就是數據之間存在着線性關系，但是在實際生活中，很多數據之間是非線性關系，雖然也可以用線性回歸擬合非線性回歸，但是效果將會很差。這個時候可以嘗試使用多項式回歸。多項式回歸中，加入了特征的更高次方（例如平方項或立方項），也相當於增加了模型的自由度，用來捕獲數據中非線性的變化。添加高階項的時候，也增加了模型的復雜度。隨着模型復雜度的升高，模型的容量以及擬合數據的能力增加，可以進一步降低訓練誤差，但導致過擬合的風險也隨之增加。
圖1 線性回歸與多項式回歸
如圖1所示，左圖數據呈現出線性關系，此時用線性回歸可以得到較好的擬合效果，如圖中的藍色直線。右圖數據則呈現出非線性的關系，顯然直線不能得到好的擬合效果，可以通過多項式回歸模型得到較好的擬合結果。
因為圖中的數據是人工生成的，所以右圖中的綠色曲線模型參數是已知的，其多項式方程為\(y=0.5x^2+x+1\)。
1. 多項式回歸的一般形式
在多項式回歸中，最重要的參數是最高次方的次數。假設最高次方的次數是\(n\)，且樣本只有一個特征時，其多項式回歸方程表示為：
如果令\(x_0=1\)，在多樣本的情況下，可以寫成向量化的形式：
其中\(X\)是大小為\(m\times (n+1)\)的矩陣，\(\theta\)是大小為\((n+1)\times 1\)的矩陣。在這里雖然只有一個特征\(x\)以及\(x\)的不同次方，但是也可以將\(x\)的高次方當做一個新特征。與多元回歸分析唯一不同的是，這些特征之間是高度相關的，而不是通常要求的那樣是相互對立的。
如果樣本不止一個特征，以只有兩個特征的2次方多項式的回歸模型為例：
其他多特征多次方的方程以此類推即可。
在這里有個問題，如果假設中出現了高階項，那么這個模型還是線性模型嗎？此時看待問題的角度不同，得到的結果也不同。如果把上面的假設看成是特征\(x\)的方程，那么該方程就是非線性方程；如果看成是參數\(\theta\)的方程，那么\(x\)的高階項都可以看做是對應\(\theta\)的參數，那么該方程就是線性方程。很明顯，在線性回歸中采用了后一種解釋方式。因此多項式回歸仍然是參數的線性模型。
回到線性模型，\(h_θ(\pmb x)=θ_0+θ_1x_1+...+θ_nx_n\), 如果我們令式\((3)\)中的\(x_0=1,x_1=x_1,x_2=x_2,x_3=x^2_1,x_4=x_2^2,x_5=x_1x_2\) ,這樣我們就得到了下式：
可以發現，我們又重新回到了線性回歸，這是一個五元線性回歸，可以用線性回歸的方法來完成算法。對於每個二元樣本特征\((x1,x2)\),我們得到一個五元樣本特征\((1,x_1,x_2,x^2_1,x_2^2,x_1x_2)\)，通過這個改進的五元樣本特征，我們把非線性回歸的函數變回了線性回歸。
2. 多項式回歸示例
下面主要使用了numpy、matplotlib和scikit-learn，如下：
# python 3.6
# sklearn 0.20.3
import numpy as np
import matplotlib.pyplot as plt
import matplotlib as mpl
from sklearn.linear_model import LinearRegression
from sklearn.preprocessing import PolynomialFeatures
from sklearn.metrics import mean_squared_error
from sklearn.pipeline import Pipeline
接下來生成只有一個特征，共10個樣本的訓練數據，使用\(y=0.5x^2+x+2\)並加入一些隨機誤差生成對應的輸出值。
rng = np.random.RandomState(1)
# 生成樣本數據
x = np.linspace(-3, 7, 10)
y = 0.5 * x**2 + x + 2 + rng.normal(0, 1.2, 10)
# 轉換成矩陣的形式
X = x[:,np.newaxis]
y = y[:,np.newaxis]
# 將樣本點顯示在坐標軸上，如圖2所示
plt.figure(figsize=(8, 6))
plt.scatter(X, y, color="r")
plt.xlabel("X")
plt.ylabel("y")
# plt.show()
圖2 生成的原始數據
2.1 使用直線方程擬合
下面先使用直線方程擬合生成的數據：
lin_reg = LinearRegression()
lin_reg.fit(X, y)
print(lin_reg.intercept_, lin_reg.coef_)
# [5.05753609] [[2.95924372]]
# 繪制出擬合數據得到的直線，如圖3所示
X_plot = np.linspace(-3, 7, 200).reshape(-1, 1)
y_plot = lin_reg.coef_*X_plot + lin_reg.intercept_
plt.figure(figsize=(8, 6))
plt.plot(X_plot, y_plot, color='b', linewidth=2)
plt.scatter(X, y, color='r')
plt.xlabel('X')
plt.ylabel('y')
# plt.show()
圖3 直線擬合效果
可以使用函數"mean_squared_error"來計算模型在訓練數據上的均方誤差(MSE):
print("線性回歸模型的MSE：{}".format(mean_squared_error(y, lin_reg.predict(X))))
# 線性回歸模型的MSE：25.902093991686804
2.2 使用多項式回歸
為了擬合2次方程，需要有特征\(x^2\)的數據，這里可以使用函數"PolynomialFeatures"來獲得：
# 參數degree設定2次方項
poly_features = PolynomialFeatures(degree=2, include_bias=False)
X_poly = poly_features.fit_transform(X)
# print(X_poly)
# 訓練模型
lin_reg = LinearRegression()
lin_reg.fit(X_poly, y)
print(lin_reg.intercept_, lin_reg.coef_)
# [1.6761218] [[0.77246083 0.54669572]]
# 將擬合得到的曲線繪制出來，如圖4所示
X_plot = np.linspace(-3, 7, 200).reshape(-1, 1)
X_plot_poly = poly_features.fit_transform(X_plot)
y_plot = np.dot(X_plot_poly, lin_reg.coef_.T) + lin_reg.intercept_
plt.figure(figsize=(8, 6))
plt.plot(X_plot, y_plot, color='b', linewidth=2)
plt.scatter(X, y, color='r')
plt.xlabel('X')
plt.ylabel('y')
# plt.show()
圖4 多項式擬合數據
print("多項式回歸模型的MSE：{:.3f}".format(mean_squared_error(y, lin_reg.predict(X_poly))))
多項式回歸模型的MSE：1.850
從上述代碼中，我們得到了多項式回歸模型的參數，即多項式方程為\(h(x)=0.55x^2+0.77x+1.68\)，而我們生成數據使用的方程為\(y=0.5x^2+x+2\)。可以看出兩個方程之間還是比較相近的。
如上所示，利用多項式回歸，損失函數MSE的值下降到了1.85。通過觀察代碼，可以發現訓練多項式方程與直線方程唯一的差別是輸入的訓練集\(X\)的差別。在訓練直線方程時直接輸入了\(X\)的值，在訓練多項式方程的時候，還添加了我們計算出來的\(x^2\)這個“新特征”的值（由於\(x^2\)完全是由\(x\)的值確定的，因此嚴格意義上來講此時該模型只有一個特征\(x\)）。
此時有個非常有趣的問題：假如一開始得到的數據就是上面代碼中"X_poly"的樣子，且不知道\(x_1\)與\(x_2\)之間的關系。此時相當於我們有10個樣本，每個樣本具有\(x_1,x_2\)兩個不同的特征。這時假設函數為：
直接按照二元線性回歸方程來訓練，也可以得到上面同樣的結果（\(\theta\)的值）。如果在相同情況下，收集到了新的數據，可以直接帶入上面的方程進行預測。唯一不同的是，我們不知道\(x_2=x^2_1\)這個隱含在數據內部的關系，所有也就無法畫出圖4中的這條曲線。一旦了解到了這兩個特征之間的關系，數據的維度就從3維下降到了2維（包含截距項\(\theta_0\)）。
2.3 多項式階數與過擬合
在上面實現多項式回歸的過程中，通過引入高階項\(x^2\)，訓練誤差從25.90下降到了1.85，減小了非常多。那么訓練誤差是否還有進一步下降的空間呢？答案是肯定的，通過繼續增加更高階的項，訓練誤差可以進一步降低。通過嘗試，當最高階項為\(x^9\)時，訓練誤差幾乎等於0了。
下面是測試不同階數(degree)產生的擬合效果：
import numpy as np
import matplotlib.pyplot as plt
import matplotlib as mpl
from sklearn.linear_model import LinearRegression
from sklearn.preprocessing import PolynomialFeatures
from sklearn.metrics import mean_squared_error
from sklearn.pipeline import Pipeline
mpl.rcParams['font.sans-serif'] = ['SimHei']
mpl.rcParams['axes.unicode_minus'] = False
np.set_printoptions(suppress=True)
rng = np.random.RandomState(1)
N = 10
x = np.linspace(-3, 7, N)
x = np.sort(x)
y = 0.5 * x**2 + x + 2 + rng.normal(0, 1.2, N)
x.shape = -1, 1
y.shape = -1, 1
model = Pipeline([
('poly', PolynomialFeatures()),
('linear', LinearRegression())])
plt.figure(figsize=(8, 6), facecolor='w')
# 多項式的階數
d_pool = np.arange(1, N, 1)
m = d_pool.size
# 配置顏色
clrs = []
rdbu = mpl.cm.get_cmap("RdYlBu", m)
for c in np.linspace(0, 1, m):
clrs.append(rdbu(c))
plt.scatter(x, y, color="r")
for i, d in enumerate(d_pool):
model.set_params(poly__degree=d)
model.fit(x, y)
# lin = model.get_params('linear')['linear']
x_hat = np.linspace(x.min(), x.max(), num=100)
x_hat.shape = -1, 1
y_hat = model.predict(x_hat)
s = mean_squared_error(y, model.predict(x))
plt.plot(x_hat, y_hat, color=clrs[i], lw=2, label=('%d階，train mse=%.3f' % (d, s)))
plt.legend(loc='upper left')
plt.grid(True)
plt.xlabel('X', fontsize=14)
plt.ylabel('Y', fontsize=14)
plt.tight_layout(1, rect=(0, 0, 1, 0.95))
plt.show()
此時，我們繪制出不同多項式階數擬合效果的圖像，如圖5所示：
圖5 不同階數擬合數據
由圖5可以看到，當多項式階數逐漸增加時，訓練誤差逐漸減低。當階數為9時，函數圖像幾乎穿過了每一個樣本點，所有的訓練樣本都落在了擬合的曲線上，訓練誤差接近於0。 可以說是近乎完美的模型了。但是，這樣的曲線與我們最開始數據的來源（一個二次方程加上一些隨機誤差）差異非常大。如果從相同來源再取一些樣本點，使用該模型預測會出現非常大的誤差。類似這種訓練誤差非常小，但是新數據點的測試誤差非常大的情況，就叫做模型的過擬合。過擬合出現時，表示模型過於復雜，過多考慮了當前樣本的特殊情況以及噪音（模型學習到了當前訓練樣本非全局的特性），使得模型的泛化能力下降。
出現過擬合一般有以下幾種解決方式：
- 降低模型復雜度，例如減小上面例子中的degree；
- 降維，減小特征的數量；
- 增加訓練樣本，在現實機器學習任務中很難辦到；
- 添加正則化項.
防止模型過擬合是機器學習領域里最重要的問題之一。鑒於該問題的普遍性和重要性，在滿足要求的情況下，能選擇簡單模型時應該盡量選擇簡單的模型。
參考來源：
1）https://www.cnblogs.com/Belter/p/8530222.html
2）https://blog.csdn.net/qq_25560849/article/details/80543180
3）https://www.cnblogs.com/pinard/p/6004041.html