scikit-learn線性迴歸，多元迴歸，多項式迴歸的實現
匹薩的直徑與價格的資料
%matplotlib inline import matplotlib.pyplot as plt def runplt(): plt.figure() plt.title(u'diameter-cost curver') plt.xlabel(u'diameter') plt.ylabel(u'cost') plt.axis([0,25,25]) plt.grid(True) return plt plt = runplt() X = [[6],[8],[10],[14],[18]] y = [[7],[9],[13],[17.5],[18]] plt.plot(X,y,'k.') plt.show()
訓練模型
from sklearn.linear_model import LinearRegression
import numpy as np
# 建立並擬合模型
model = LinearRegression()
model.fit(X,y)
print('預測一張12英寸匹薩價格：$%.2f' % model.predict(np.array([12]).reshape(-1,1))[0])
預測一張12英寸匹薩價格：$13.68
一元線性迴歸假設解釋變數和響應變數之間存線上性關係；這個線性模型所構成的空間是一個超平面（hyperplane）。
超平面是n維歐氏空間中餘維度等於一的線性子空間，如平面中的直線、空間中的平面等，總比包含它的空間少一維。
在一元線性迴歸中，一個維度是響應變數，另一個維度是解釋變數，總共兩維。因此，其超平面只有一維，就是一條線。
上述程式碼中sklearn.linear_model.LinearRegression類是一個估計器（estimator）。估計器依據觀測值來預測結果。在scikit-learn裡面，所有的估計器都帶有:
- fit()
- predict()
fit()用來分析模型引數，predict()是通過fit()算出的模型引數構成的模型，對解釋變數進行預測獲得的值。
因為所有的估計器都有這兩種方法，所有scikit-learn很容易實驗不同的模型。
一元線性迴歸模型：
y=α+βx
一元線性迴歸擬合模型的引數估計常用方法是:
- 線性最小二乘法（linear least squares）
首先，我們定義出擬合成本函式，然後對引數進行數理統計。
plt = runplt() plt.plot(X,'k.') X2 = [[0],[25]] model = LinearRegression() model.fit(X,y) y2 = model.predict(X2) plt.plot(X,'k.') plt.plot(X2,y2,'g-') plt.show()
plt = runplt() plt.plot(X,'k.') y3 = [14.25,14.25,14.25] y4 = y2 * 0.5 + 5 model.fit(X[1:-1],y[1:-1]) y5 = model.predict(X2) plt.plot(X,'g-.') plt.plot(X2,y3,'r-.') plt.plot(X2,y4,'y-.') plt.plot(X2,y5,'o-') plt.show()
成本函式（cost function）也叫損失函式（loss function），用來定義模型與觀測值的誤差。模型預測的價格與訓練集資料的差異稱為殘差（residuals）或訓練誤差（training errors）。後面我們會用模型計算測試集，那時模型預測的價格與測試集資料的差異稱為預測誤差（prediction errors）或訓練誤差（test errors）。
模型的殘差是訓練樣本點與線性迴歸模型的縱向距離，如下圖所示：
plt = runplt() plt.plot(X,'g-') # 殘差預測值 yr = model.predict(X) for idx,x in enumerate(X): plt.plot([x,x],[y[idx],yr[idx]],'r-') plt.show()
我們可以通過殘差之和最小化實現最佳擬合，也就是說模型預測的值與訓練集的資料最接近就是最佳擬合。對模型的擬合度進行評估的函式稱為殘差平方和（residual sum of squares）成本函式。就是讓所有訓練資料與模型的殘差的平方之和最小化，如下所示：
其中，
import numpy as np
print('殘差平方和: %.2f' % np.mean((model.predict(X) - y) ** 2))
殘差平方和: 1.75
解一元線性迴歸的最小二乘法
通過成本函式最小化獲得引數，我們先求相關係數 ββ 。按照頻率論的觀點，我們首先需要計算 xx 的方差和 xx 與 yy 的協方差。
方差是用來衡量樣本分散程度的。如果樣本全部相等，那麼方差為0。方差越小，表示樣本越集中，反正則樣本越分散。方差計算公式如下：
Numpy裡面有var方法可以直接計算方差，ddof引數是貝塞爾(無偏估計)校正係數（Bessel's correction），設定為1，可得樣本方差無偏估計量。
print(np.var([6,8,10,14,18],ddof=1))
23.2
協方差表示兩個變數的總體的變化趨勢。如果兩個變數的變化趨勢一致，也就是說如果其中一個大於自身的期望值，另外一個也大於自身的期望值，那麼兩個變數之間的協方差就是正值。 如果兩個變數的變化趨勢相反，即其中一個大於自身的期望值，另外一個卻小於自身的期望值，那麼兩個變數之間的協方差就是負值。如果兩個變數不相關，則協方差為0，變數線性無關不表示一定沒有其他相關性。協方差公式如下：
其中，
import numpy as np print(np.cov([6,[7,9,13,17.5,18])[0][1])
22.65
現在有了方差和協方差，就可以計算相關係統
算出
將前面的資料帶入公式就可以求出
模型評估
前面我們用學習演算法對訓練集進行估計，得出了模型的引數。有些度量方法可以用來評估預測效果，我們用R方（r-squared）評估匹薩價格預測的效果。R方也叫確定係數（coefficient of determination），表示模型對現實資料擬合的程度。計算R方的方法有幾種。一元線性迴歸中R方等於皮爾遜積矩相關係數（Pearson product moment correlation coefficient或Pearson's r）的平方。種方法計算的R方一定介於0～1之間的正數。其他計算方法，包括scikit-learn中的方法，不是用皮爾遜積矩相關係數的平方計算的，因此當模型擬合效果很差的時候R方會是負值。下面我們用scikit-learn方法來計算R方。
R方是0.6620說明測試集裡面過半數的價格都可以通過模型解釋。現在，讓我們用scikit-learn來驗證一下。LinearRegression的score方法可以計算R方：
# 測試集 X_test = [[8],[11],[16],[12]] y_test = [[11],[8.5],[15],[18],[11]] model = LinearRegression() model.fit(X,y) model.score(X_test,y_test)
0.66200528638545164
多元迴歸
from sklearn.linear_model import LinearRegression
X = [[6,2],[8,1],[10,0],[14,[18,0]]
y = [[7],[18]]
model = LinearRegression()
model.fit(X,y)
X_test = [[8,[9,[11,[16,[12,0]]
y_test = [[11],[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))
Predicted: [ 10.06250019],Target: [11]
Predicted: [ 10.28125019],Target: [8.5]
Predicted: [ 13.09375019],Target: [15]
Predicted: [ 18.14583353],Target: [18]
Predicted: [ 13.31250019],Target: [11]
R-squared: 0.77
多項式迴歸
上例中，我們假設解釋變數和響應變數的關係是線性的。真實情況未必如此。下面我們用多項式迴歸，一種特殊的多元線性迴歸方法，增加了指數項。現實世界中的曲線關係都是通過增加多項式實現的，其實現方式和多元線性迴歸類似。本例還用一個解釋變數，匹薩直徑。讓我們用下面的資料對兩種模型做個比較：
import numpy as np
from sklearn.linear_model import LinearRegression
from sklearn.preprocessing import PolynomialFeatures
X_train = [[6],[18]]
y_train = [[7],[18]]
X_test = [[6],[16]]
y_test = [[8],[12],[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 = runplt()
plt.plot(X_train,y_train,'k.')
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),'r-')
plt.show()
print(X_train)
print(X_train_quadratic)
print(X_test)
print(X_test_quadratic)
print('1 r-squared',regressor.score(X_test,y_test))
print('2 r-squared',regressor_quadratic.score(X_test_quadratic,y_test))
[[6],[18]]
[[ 1. 6. 36.]
[ 1. 8. 64.]
[ 1. 10. 100.]
[ 1. 14. 196.]
[ 1. 18. 324.]]
[[6],[16]]
[[ 1. 6. 36.]
[ 1. 8. 64.]
[ 1. 11. 121.]
[ 1. 16. 256.]]
('1 r-squared',0.80972683246686095)
('2 r-squared',0.86754436563450732)
plt = runplt()
plt.plot(X_train,'k.')
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,'r-')
cubic_featurizer = PolynomialFeatures(degree=3)
X_train_cubic = cubic_featurizer.fit_transform(X_train)
X_test_cubic = cubic_featurizer.transform(X_test)
regressor_cubic = LinearRegression()
regressor_cubic.fit(X_train_cubic,y_train)
xx_cubic = cubic_featurizer.transform(xx.reshape(xx.shape[0],regressor_cubic.predict(xx_cubic))
plt.show()
print(X_train_cubic)
print(X_test_cubic)
print('2 r-squared',y_test))
print('3 r-squared',regressor_cubic.score(X_test_cubic,y_test))
[[ 1.00000000e+00 6.00000000e+00 3.60000000e+01 2.16000000e+02]
[ 1.00000000e+00 8.00000000e+00 6.40000000e+01 5.12000000e+02]
[ 1.00000000e+00 1.00000000e+01 1.00000000e+02 1.00000000e+03]
[ 1.00000000e+00 1.40000000e+01 1.96000000e+02 2.74400000e+03]
[ 1.00000000e+00 1.80000000e+01 3.24000000e+02 5.83200000e+03]]
[[ 1.00000000e+00 6.00000000e+00 3.60000000e+01 2.16000000e+02]
[ 1.00000000e+00 8.00000000e+00 6.40000000e+01 5.12000000e+02]
[ 1.00000000e+00 1.10000000e+01 1.21000000e+02 1.33100000e+03]
[ 1.00000000e+00 1.60000000e+01 2.56000000e+02 4.09600000e+03]]
('2 r-squared',0.86754436563450732)
('3 r-squared',0.83569241560369567)
plt = runplt()
plt.plot(X_train,'r-')
seventh_featurizer = PolynomialFeatures(degree=7)
X_train_seventh = seventh_featurizer.fit_transform(X_train)
X_test_seventh = seventh_featurizer.transform(X_test)
regressor_seventh = LinearRegression()
regressor_seventh.fit(X_train_seventh,y_train)
xx_seventh = seventh_featurizer.transform(xx.reshape(xx.shape[0],regressor_seventh.predict(xx_seventh))
plt.show()
print('2 r-squared',y_test))
print('7 r-squared',regressor_seventh.score(X_test_seventh,y_test))
('2 r-squared',0.86754436563450732)
('7 r-squared',0.49198460568655)
可以看出，七次擬合的R方值更低，雖然其圖形基本經過了所有的點。可以認為這是擬合過度（over-fitting）的情況。這種模型並沒有從輸入和輸出中推匯出一般的規律，而是記憶訓練集的結果，這樣在測試集的測試效果就不好了。
正則化
LASSO方法會產生稀疏引數，大多數相關係數會變成0，模型只會保留一小部分特徵。而嶺迴歸還是會保留大多數儘可能小的相關係數。當兩個變數相關時，LASSO方法會讓其中一個變數的相關係數會變成0，而嶺迴歸是將兩個係數同時縮小。
import numpy as np
from sklearn.datasets import load_boston
from sklearn.linear_model import SGDRegressor
from sklearn.cross_validation import cross_val_score
from sklearn.preprocessing import StandardScaler
from sklearn.cross_validation import train_test_split
data = load_boston()
X_train,X_test,y_test = train_test_split(data.data,data.target)
X_scaler = StandardScaler()
y_scaler = StandardScaler()
X_train = X_scaler.fit_transform(X_train)
y_train = y_scaler.fit_transform(y_train.reshape(-1,1))
X_test = X_scaler.transform(X_test)
y_test = y_scaler.transform(y_test.reshape(-1,1))
regressor = SGDRegressor(loss='squared_loss',penalty="l1")
scores = cross_val_score(regressor,X_train,y_train.reshape(-1,1),cv=5)
print('cv R',scores)
print('mean of cv R',np.mean(scores))
regressor.fit_transform(X_train,y_train)
print('Test set R',y_test))
('cv R',array([ 0.74761441,0.62036841,0.6851797,0.63347999,0.79476346]))
('mean of cv R',0.69628119572104885)
('Test set R',0.75084948718041566)
以上就是本文的全部內容，希望對大家的學習有所幫助，也希望大家多多支援我們。