本文將用一個例子來講述怎么用scikit-learn和pandas來學習Ridge回歸。
1. Ridge回歸的損失函數
在我的另外一遍講線性回歸的文章中，對Ridge回歸做了一些介紹，以及什么時候適合用 Ridge回歸。如果對什么是Ridge回歸還完全不清楚的建議閱讀我這篇文章。
Ridge回歸的損失函數表達形式是：
\(J(\mathbf\theta) = \frac{1}{2}(\mathbf{X\theta} - \mathbf{Y})^T(\mathbf{X\theta} - \mathbf{Y}) + \frac{1}{2}\alpha||\theta||_2^2\)
其中\(\alpha\)為常數系數，需要進行調優。\(||\theta||_2\)為L2范數。
算法需要解決的就是在找到一個合適的超參數\(\alpha\)情況下，求出使\(J(\mathbf\theta)\)最小的\(\theta\)。一般可以用梯度下降法和最小二乘法來解決這個問題。scikit-learn用的是最小二乘法。
2. 數據獲取與預處理
這里我們仍然用UCI大學公開的機器學習數據來跑Ridge回歸。
數據的介紹在這： http://archive.ics.uci.edu/ml/datasets/Combined+Cycle+Power+Plant
數據的下載地址在這： http://archive.ics.uci.edu/ml/machine-learning-databases/00294/
完整的代碼見我的github: https://github.com/ljpzzz/machinelearning/blob/master/classic-machine-learning/ridge_regression_1.ipynb
里面是一個循環發電場的數據，共有9568個樣本數據，每個數據有5列，分別是:AT（溫度）, V（壓力）, AP（濕度）, RH（壓強）, PE（輸出電力)。我們不用糾結於每項具體的意思。
我們的問題是得到一個線性的關系，對應PE是樣本輸出，而AT/V/AP/RH這4個是樣本特征， 機器學習的目的就是通過調節超參數\(\alpha\)得到一個線性回歸模型，即:
\(PE = \theta_0 + \theta_1*AT + \theta_2*V + \theta_3*AP + \theta_4*RH\)
使損失函數\(J(\mathbf\theta)\)最小。而需要學習的，就是\(\theta_0, \theta_1, \theta_2, \theta_3, \theta_4\)這5個參數。
下載后的數據可以發現是一個壓縮文件，解壓后可以看到里面有一個xlsx文件，我們先用excel把它打開，接着“另存為“”csv格式，保存下來，后面我們就用這個csv來運行Ridge回歸。
這組數據並不一定適合用Ridge回歸模型，實際上這組數據是高度線性的，使用正則化的Ridge回歸僅僅只是為了講解方便。
3. 數據讀取與訓練集測試集划分
我們先打開ipython notebook,新建一個notebook。當然也可以直接在python的交互式命令行里面輸入，不過還是推薦用notebook。下面的例子和輸出我都是在notebook里面跑的。
先把要導入的庫聲明了：
import matplotlib.pyplot as plt %matplotlib inline import numpy as np import pandas as pd from sklearn import datasets, linear_model
接着用pandas讀取數據：
# read_csv里面的參數是csv在你電腦上的路徑，此處csv文件放在notebook運行目錄下面的CCPP目錄里 data = pd.read_csv('.\CCPP\ccpp.csv')
我們用AT， V，AP和RH這4個列作為樣本特征。用PE作為樣本輸出：
X = data[['AT', 'V', 'AP', 'RH']] y = data[['PE']]
接着把數據集划分為訓練集和測試集：
from sklearn.cross_validation import train_test_split X_train, X_test, y_train, y_test = train_test_split(X, y, random_state=1)
4. 用scikit-learn運行Ridge回歸
要運行Ridge回歸，我們必須要指定超參數\(\alpha\)。你也許會問：“我也不知道超參數是多少啊？” 我也不知道，那么我們隨機指定一個(比如1)，后面我們會講到用交叉驗證從多個輸入超參數\(\alpha\)中快速選擇最優超參數的辦法。
from sklearn.linear_model import Ridge ridge = Ridge(alpha=1) ridge.fit(X_train, y_train)
訓練完了，可以看看模型參數是多少:
print ridge.coef_ print ridge.intercept_
輸出結果如下：
也就是說我們得到的模型是：
\(PE = 447.05552892 - 1.97373209*AT - 0.2323016*V + 0.06935852*AP - 0.15806479*RH\)
但是這樣還沒有完？為什么呢，因為我們假設了超參數\(\alpha\)為1， 實際上我們並不知道超參數\(\alpha\)取多少最好，實際研究是需要在多組自選的\(\alpha\)中選擇一個最優的。
那么我們是不是要把上面這段程序在N種\(\alpha\)的值情況下，跑N遍，然后再比較結果的優劣程度呢？ 可以這么做，但是scikit-learn提供了另外一個交叉驗證選擇最優\(\alpha\)的API，下面我們就用這個API來選擇\(\alpha\)。
5. 用scikit-learn選擇Ridge回歸超參數\(\alpha\)
這里我們假設我們想在這10個\(\alpha\)值中選擇一個最優的值。代碼如下：
from sklearn.linear_model import RidgeCV ridgecv = RidgeCV(alphas=[0.01, 0.1, 0.5, 1, 3, 5, 7, 10, 20, 100]) ridgecv.fit(X_train, y_train) ridgecv.alpha_
輸出結果為：7.0，說明在我們給定的這組超參數中， 7是最優的\(\alpha\)值。
6. 用scikit-learn研究超參數\(\alpha\)和回歸系數\(\theta\)的關系
通過Ridge回歸的損失函數表達式可以看到，\(\alpha\)越大，那么正則項懲罰的就越厲害，得到回歸系數\(\theta\)就越小，最終趨近與0。而如果\(\alpha\)越小，即正則化項越小，那么回歸系數\(\theta\)就越來越接近於普通的線性回歸系數。
這里我們用scikit-learn來研究這種Ridge回歸的變化，例子參考了scikit-learn的官網例子。我們單獨啟動一個notebook或者python shell來運行這個例子。
完整的代碼見我的github: https://github.com/ljpzzz/machinelearning/blob/master/classic-machine-learning/ridge_regression.ipynb
首先還是加載類庫：
import numpy as np import matplotlib.pyplot as plt from sklearn import linear_model %matplotlib inline
接着我們自己生成一個10x10的矩陣X，表示一組有10個樣本，每個樣本有10個特征的數據。生成一個10x1的向量y代表樣本輸出。
# X is a 10x10 matrix X = 1. / (np.arange(1, 11) + np.arange(0, 10)[:, np.newaxis]) # y is a 10 x 1 vector y = np.ones(10)
這樣我們的數據有了，接着就是准備超參數\(\alpha\)了。我們准備了200個超參數，來分別跑 Ridge回歸。准備這么多的目的是為了后面畫圖看\(\alpha\)和\(\theta\)的關系
n_alphas = 200 # alphas count is 200, 都在10的-10次方和10的-2次方之間 alphas = np.logspace(-10, -2, n_alphas)
有了這200個超參數\(\alpha\)，我們做200次循環，分別求出各個超參數對應的\(\theta\)(10個維度)，存起來后面畫圖用。
clf = linear_model.Ridge(fit_intercept=False) coefs = [] # 循環200次 for a in alphas: #設置本次循環的超參數 clf.set_params(alpha=a) #針對每個alpha做ridge回歸 clf.fit(X, y) # 把每一個超參數alpha對應的theta存下來 coefs.append(clf.coef_)
好了，有了200個超參數\(\alpha\)，以及對應的\(\theta\)，我們可以畫圖了。我們的圖是以\(\alpha\)為x軸，\(\theta\)的10個維度為y軸畫的。代碼如下：
ax = plt.gca() ax.plot(alphas, coefs) #將alpha的值取對數便於畫圖 ax.set_xscale('log') #翻轉x軸的大小方向，讓alpha從大到小顯示 ax.set_xlim(ax.get_xlim()[::-1]) plt.xlabel('alpha') plt.ylabel('weights') plt.title('Ridge coefficients as a function of the regularization') plt.axis('tight') plt.show()
最后得到的圖如下：
從圖上也可以看出，當\(\alpha\)比較大，接近於\(10^{-2}\)的時候，\(\theta\)的10個維度都趨於0。而當\(\alpha\)比較小，接近於\(10^{-10}\)的時候，\(\theta\)的10個維度都趨於線性回歸的回歸系數。
（歡迎轉載，轉載請注明出處。歡迎溝通交流： liujianping-ok@163.com）