[Python] ML WITH SCIKIT Learn (Model Selection)

年初網路上有一個基礎 scikit learning 的教學, 相當的完整以及實用

花了一點時間看完發覺應該讓更多人學會使用他, 所以打算來寫個中文版本的說明

這篇只是用我的話重新講述 Model Selection 部份, 剩餘的有時間再寫

Video: https://www.youtube.com/watch?v=iFkRt3BCctg

官網: https://us.pycon.org/2013/community/tutorials/23/

影片相關資料: https://github.com/ogrisel/parallel_ml_tutorial.git  (需要用 git 下載)

下面程式碼一律都是 Olivier Grisel 所提供在 github 的, 我只是會有小幅度的修改跟利用

想要看完整版還是乖乖的把 video 看完吧, 要花約 2.5 hr (Model selection只佔 1/4)

Model Selection

首先如果想要拿 python 做一些 data analysis 或者是 machine learning 的應用

最好手上先有一些工具使用會比較方便

首先是常見要先準備好的 Tool :

numpy, scipy(皆為矩陣運算或者數學運算較為方便)
(http://www.scipy.org/)

pylab (其實他包含有Matplotlib 繪圖, 這邊是拿來畫分布圖或者是直條圖之類)
(http://wiki.scipy.org/PyLab)

sklearn (machine learning library, 基本上以使用上最為簡單的 ML library)
(http://scikit-learn.org/stable/)

上述的工具最好都要先裝好, 後續會比較方便

由於我是用 Ubuntu, 所以安裝不外乎 apt-get 或者是安裝好 pip, 用 pip install 安裝 python library

基本上整個流程可以濃縮程成幾個步驟

首先 import 下面的 library

import numpy as np
import pylab as pl

Prepare data

首先針對你要處理的 data, 要先有一定的認知

假設我今天用 scikit 的範例 data - Digits
(這篇會用這個 dataset 做一些範例做說明, 而這是一個辨識數字的 Database)

from sklearn.datasets import load_digits
digits = load_digits()
X, y = digits.data, digits.target

print("data shape: %r, target shape: %r" % (X.shape, y.shape))
print("classes: %r" % list(np.unique(y)))

n_samples, n_features = X.shape

print("n_samples=%d" % n_samples)
print("n_features=%d" % n_features)

因為內建的 dataset scikit都幫你整理好了, 所以可以直接使用一些變數

像是 digits.data 就會是真實的 data  資料, digits.target 就是要被 training 的 classes

那因為 digits.data 實際上是 numpy 的 array, 所以可以用 shape 指令找出他的維度

印出結果

data shape: (1797, 64), target shape: (1797,)
classes: [0, 1, 2, 3, 4, 5, 6, 7, 8, 9]
n_samples=1797
n_features=64

事實上, 就算今天用的不是他的 dataset, 也應該要先找出你要用的 data 相關的資料

舉凡是檔案的維度, 這邊來說也就是 [筆數] X [feature 數量]

對應的 class  數量正確與否,  class 有多少種類別, 檔案建立的時間

甚至是是否有 missing value, 每個 feature 的類別都應該有紀錄(數字, 文字...etc)

這樣對後續的實驗才會有比較完整的概要

Draw - Randomized PCA

接下來是針對 data 做一些 PCA 的動作, 甚麼? 這甚麼鬼? 為甚麼跳那麼快?

其實因為當你拿到一堆 multi-class 的 data, 可能很難馬上了解有甚麼特性可以用

這時候可以看 data 的狀況, 因為他用的是有 10 個 class 的 data

所以可以先用 PCA 幫忙篩選出重要的 feature 然後繪圖成 2D 來輔助

(為甚麼 PCA 可以篩選出重要的 feature? 請參看台大張智星教授的PCA教學影片解釋
 http://www.youtube.com/watch?v=pAuK4bayhdw )

流程如下

首先, 先利用 sklearn 的 RandomizedPCA 將 data 降維到剩 2 維

(簡單來說 PCA 是會利用數學上的 eigenvector, eigenvalue 屬性找出 data 中比較重要的維度出來並且依序排列, 所以這邊會有找出最重要的兩維的含意)

from sklearn.decomposition import RandomizedPCA

pca = RandomizedPCA(n_components=2)
X_pca = pca.fit_transform(X)

print X_pca.shape  # print(1797, 2), 2 dimension left

然後就是畫畫的時間啦!! (一樣, 用他範例的 code 呈現)

from itertools import cycle

colors = ['b', 'g', 'r', 'c', 'm', 'y', 'k']
markers = ['+', 'o', '^', 'v', '<', '>', 'D', 'h', 's']
print y
for i, c, m in zip(np.unique(y), cycle(colors), cycle(markers)):
    pl.scatter(X_pca[y == i, 0], X_pca[y == i, 1], 
        c=c, marker=m, label=i, alpha=0.5)
    
pl.legend(loc='best')
pl.show()

語法不熟的可以參照 http://matplotlib.org/ 網站

結果會顯示出下圖

看起來五顏六色齁, 不過要是仔細看, 每個顏色代表不同的 class

會發現如果只看兩個重要的維度, 好像並沒有太難分割出彼此的類別, 或許或多或少有點差異

但是應該還算是可以區分出來的

當然因為他範例拿的比較剛好, 並非所有 data 都可以這麼剛好畫出這樣的結果

首先因為 class 夠多, 所以這樣呈現比較相對清楚

再來因為其特徵差異其實已經有一定程度的不同, 所以更加強了這顯示結果

我拿一個比較悽慘的例子, 下圖也是做過 RPCA 的結果

這是我自己有在分析的 data, 由於 feature 不多且彼此差異不算明顯, 只有兩個 class 類別

所以會看到, 就算是最重要的 feature, 彼此都已經互相交疊很嚴重了

這時候用這方式呈現就挺悽慘的

因為這是我自己找出來的 feature, 通常我會重新設計或者是找尋更好的 feature

Train

就以一個完全不懂 Machine Learning 的人來說, 又想要直接用 ML library 最簡單的作法

就是讀取檔案之後, 直接跑演算法 (假設檔案不大跑得動)

演算法很難 implement ?

罰你去上台大資工林軒田教授的 Machine Learning 課, 快快樂樂實作 ML 馬上就會啦!!(大誤)

其實 scikit 內建超多的 ML 演算法, 交給 scikit learn 就好啦!!

馬上就拿傳說中的 ML 演算法試試看 - SVC
(註: scikit-learn SVC 是使用 libsvm 的模組
libsvm 是台大林智仁老師實驗室所寫的http://www.csie.ntu.edu.tw/~cjlin/libsvm/)

from sklearn.svm import SVC
print SVC().fit(X, y).score(X, y) # print 1.0

在 scikit learn 裡面 SVC() 就是 support vector classification 演算法

fit() 就是 training 的意思, score(X, y) 的意思就是給定 X data, 我去預估 X 的結果

然後會跟你給定的 y 做比較, 然後給出命中率

那....一般人看到結果可能超高興的, 居然是 1.0, 1.0 就是 100%, 萬歲!!!

但是...別高興的太早, 基本上有一定程度了解 ML 得人心中一定會馬上冒出一個字眼 - Overfit

Overfit 簡單講就是, 雖然我請你學習並且解數學 100 題, 結果你死背答案所以都對了

但是當我給你類似的題目但是數字不同的時候, 你就慘了

Overfit 就有點這樣的感覺, 沒有學到真的該學的

雖然 SVC 很強沒錯, 但是也無法直接確保他真的學"對"了

不過上述都只是我的假設, 如果沒有新的題目, 其實也沒辦法說出是否真的學會

所以我們來切割一下原本的 data 吧, 假設他學 8 成, 能準確預估出剩餘2成的幾成?

切割的工具不用擔心, 不用自己寫, 交給 scikit learn 就好啦!!!

scikit learn 其中有一個 cross validation library 有提供 split 功能如下

from sklearn.cross_validation import train_test_split

X_train, X_test, y_train, y_test = train_test_split(
    X, y, test_size=0.2, random_state=0)

print("train data shape: %r, train target shape: %r"
      % (X_train.shape, y_train.shape)) # train data shape: (1437, 64), train target shape: (1437,)
print("test data shape: %r, test target shape: %r"
      % (X_test.shape, y_test.shape)) # test data shape: (360, 64), test target shape: (360,)

這邊稍稍解釋一下, 我們要給 SVC 學習的就叫做 train, 要給他學完做測試的叫 test

所以這邊參數有 0.2 意思就是 20% 為 test size, 後續的 train_test_split 會幫你處理好好

切割完之後我們就有 X_train, X_test, y_train, y_test 可以用了!!

首先跟原本的方式一樣, 直接 train, 但是只放入 train

svc = SVC(kernel='rbf').fit(X_train, y_train)
train_score = svc.score(X_train, y_train) 
print train_score # 1.0 still !!

不愧是 SVC, 就算是只有 8 成的檔案也能夠準確 100%

那接下來拿來預測切割出來的 20% test

test_score = svc.score(X_test, y_test)
print test_score # 0.325 crap....

婀....好像有點悲劇, 居然只有 0.325, 也就是只有 3 成準確率, 天阿阿~~~~~~

那....八成(雙關)是 Overfit 了, 該如何是好? 通常這時候都是參數的設定有問題

而 SVC 參數預設的值是 C = 100, gamma = 0, 通常這兩個參數會大大影響整體的效果

只要經過一連串的 cv 跟 grid 之後通常可以找出不錯的設定來使用

下列是找出的結果.....我知道跳太快, 所以下一段會開始講如何用 scikit learn 做 cv

這邊先給找出適當參數的成果

svc_2 = SVC(kernel='rbf', C=10, gamma=0.001).fit(X_train, y_train)
train_score = svc_2.score(X_train, y_train) 
print train_score # 1.0 still

test_score = svc_2.score(X_test, y_test)
print test_score # 0.9916 ~ OMG

居然有有 99.1% 的準確率!! 怎麼找的? 下面會慢慢解釋~

Cross Validation (CV)

CV 其實很單純, 實驗做一次的結果通常不夠有說服力, 那我們就做多一點, data 分割的精確一點

CV 實作很麻煩? 別擔心, 統統交給 scikit learn 就好啦 ~ 哇哈哈哈哈 (筆者瘋瘋的)

scikit learn 中可以用 ShuffleSplit 來幫你切割需求不同的 data

from sklearn.cross_validation import ShuffleSplit

cv = ShuffleSplit(n_samples, n_iter=3, test_size=0.1, random_state=0)

上述的意思就是我要將原始的 data 切成 9 : 1 而且要做成 3 份, 也就是 3-fold cv

也就是如果有 1000 筆 data, 我會做出三組不完全相同的 train: 900 test:100 這樣比例的 data

放個實際使用範例

from sklearn.cross_validation import ShuffleSplit

cv = ShuffleSplit(n_samples, n_iter=3, test_size=0.1,
    random_state=0)

for cv_index, (train, test) in enumerate(cv):
    print("# Cross Validation Iteration #%d" % cv_index)
    print("train indices: {0}...".format(train[:10]))
    print("test indices: {0}...".format(test[:10]))
    
    svc = SVC(kernel="rbf", C=1, gamma=0.001).fit(X[train], y[train])
    print("train score: {0:.3f}, test score: {1:.3f}\n".format(
        svc.score(X[train], y[train]), svc.score(X[test], y[test])))

這邊看起來很玄, 其實 cv 的結果只是幫你切割好要存取的 row 成一個 list, 然後你再去存取他

結果大致上如下:

# Cross Validation Iteration #0
train indices: [ 353    5   58 1349 1025  575 1074 1110 1745  689]...
test indices: [1081 1707  927  713  262  182  303  895  933 1266]...
train score: 0.929, test score: 0.283

# Cross Validation Iteration #1
train indices: [1336  608  977   22  526 1587 1130  569 1481  962]...
test indices: [1014  755 1633  117  181  501  948 1076   45  659]...
train score: 0.866, test score: 0.850

# Cross Validation Iteration #2
train indices: [ 451  409  911 1551  133  691 1306  111  852  825]...
test indices: [ 795  697  655  573  412  743  635  851 1466 1383]...
train score: 0.863, test score: 0.883

可以看出 test 結果有好有壞!!

那你知道, scikit 是個有點萬能的 tool, 上述這麼冗的步驟, 他當然有現成的工具可以用囉

所以其實上述的步驟可以縮減成

from sklearn.cross_validation import cross_val_score

cv = ShuffleSplit(n_samples, n_iter=10, test_size=0.1, random_state=0)

test_scores = cross_val_score(svc, X, y, cv=cv, n_jobs=2) # use default svc
print test_scores # [ 0.33888889  0.52777778  0.55        0.45        0.44444444  0.57222222
  0.52777778  0.47222222  0.44444444  0.56666667]

看起來相當方便阿, 只要用 ShuffleSplit 割好 10 份, 馬上就幫你做出 10 份的 cv (test)結果!!

那你知道的, 看 vector 的 result 有點麻煩, 而且實際上好不好還要算平均啦標準差啦....

所以再偷懶一點....用 scipy 的 stats 來計算

from scipy.stats import sem

def mean_score(scores):
    """Print the empirical mean score and standard error of the mean."""
    return ("Mean score: {0:.3f} (+/-{1:.3f})").format(np.mean(scores), sem(scores))

print(mean_score(test_scores)) # Mean score: 0.489 (+/-0.023)

相當相當方便阿, 不過這 CV 有點悽慘, 不過既然知道怎麼做 CV 之後, 就可以開始來找參數了

Grid

其實就很單純, 以 SVC 來說, 有 C 跟 gamma 兩個參數要調整, 我就將很多不同的數字帶入
(其實只要是有需要調整參數的 learner 都可以先用 grid 先找出適當的參數)

然後做大量的 CV, 找平均最好的結果, 就可以啦!!

說很簡單, 來看範例最一開始的寫法

n_gammas = 10
n_iter = 5
cv = ShuffleSplit(n_samples, n_iter=n_iter, train_size=500, test_size=500, random_state=0)

train_scores = np.zeros((n_gammas, n_iter))
test_scores = np.zeros((n_gammas, n_iter))
gammas = np.logspace(-7, -1, n_gammas)

for i, gamma in enumerate(gammas):
    for j, (train, test) in enumerate(cv):
        clf = SVC(C=10, gamma=gamma).fit(X[train], y[train])
        train_scores[i, j] = clf.score(X[train], y[train])
        test_scores[i, j] = clf.score(X[test], y[test])

看起來有點複雜, 簡單講, 前面幾行都容易理解

n_gammas = 10 是指要帶入十組不同的 gamma, n_iter = 5 檔案要切割 5 種不同分佈的 data

train_score 跟 test_score 只是先宣告好等會要放結果的變數

gammas 就有趣了, 因為前面有說要帶入不同的數字, 那我先針對 gamma 來改變

那如果稍微有了解過 libsvm 的人, 會知道要帶入 gamma 的差距通常會給上下界然後抓 log 曲線
(也就是一般帶數字是 1,2,3,4 彼此差 1 他則是從 10^-7帶到10^-1, 彼此差 log 距離)

那 C 先固定, 因為現在只有要測試不同的 gamma, 這邊預設成 10 (C=10)

然後將結果存到 train_score 跟 test_score

跑完之後把結果繪圖出來看看測試情況

也就是將相對的 gamma 值跟 cv 的 test result 匯製成表

for i in range(n_iter):
    pl.semilogx(gammas, train_scores[:, i], alpha=0.4, lw=2, c='b')
    pl.semilogx(gammas, test_scores[:, i], alpha=0.4, lw=2, c='g')
pl.ylabel("score for SVC(C=10, gamma=gamma)")
pl.xlabel("gamma")
pl.text(1e-6, 0.5, "Underfitting", fontsize=16, ha='center', va='bottom')
pl.text(1e-4, 0.5, "Good", fontsize=16, ha='center', va='bottom')
pl.text(1e-2, 0.5, "Overfitting", fontsize=16, ha='center', va='bottom')
pl.show()

這邊語法要稍微熟悉一下 matplotlib 語法, semilogx 只是在繪畫的時候會將 x 軸做 log scaling

如下圖

藍線為 training result, 綠線為 test result, 藍線高不一定代表好, 綠線高才是比較好的依據

橫軸為 gamma 值, 縱軸為 CV test 的準確度, 越高越好, 所以看到最高為 0.001 (10^-3) 左右

這樣是不是就一目了然參數的不同導致不一樣的結果呢?

所以相對的, 我們現在就可以固定 gamma 為 0.001

然後開始跑不同的 C 值, 程式碼基本上跟上面差不多, 只是固定得變數換了

n_Cs = 10
n_iter = 5
cv = ShuffleSplit(n_samples, n_iter=n_iter, train_size=500, test_size=500,
    random_state=0)

train_scores = np.zeros((n_Cs, n_iter))
test_scores = np.zeros((n_Cs, n_iter))
Cs = np.logspace(-5, 5, n_Cs)

for i, C in enumerate(Cs):
    for j, (train, test) in enumerate(cv):
        clf = SVC(C=C, gamma=1e-3).fit(X[train], y[train])
        train_scores[i, j] = clf.score(X[train], y[train])
        test_scores[i, j] = clf.score(X[test], y[test])

那一樣, 把結果繪製成表

for i in range(n_iter):
    pl.semilogx(Cs, train_scores[:, i], alpha=0.4, lw=2, c='b')
    pl.semilogx(Cs, test_scores[:, i], alpha=0.4, lw=2, c='g')
pl.ylabel("score for SVC(C=C, gamma=1e-3)")
pl.xlabel("C")
pl.text(1e-3, 0.5, "Underfitting", fontsize=16, ha='center', va='bottom')
pl.text(1e3, 0.5, "Few Overfitting", fontsize=16, ha='center', va='bottom')
pl.show()

如下圖

橫軸為 C 值, 縱軸為 CV test 的準確度, 越高越好, 所以看到最高 C 大約為 10 ~ 100 的時候

所以這樣就找出適當的 C & gamma 值了

那...萬能的 scikit learn 阿, 需要做這麼多好麻煩的事情, 來點更好用的工具吧!!

沒錯, 他有提供整套的服務!! (我這篇根本已經變成幫 scikit 打廣告了XDDD)

scikit learn 有著有著完整個 grid 方法 - GridSearchCV

from sklearn.grid_search import GridSearchCV
from pprint import pprint
svc_params = {
    'C': np.logspace(-1, 2, 4),
    'gamma': np.logspace(-4, 0, 5),
}

n_subsamples = 500
X_small_train, y_small_train = X_train[:n_subsamples], y_train[:n_subsamples]

gs_svc = GridSearchCV(SVC(), svc_params, cv=10, n_jobs=-1)

gs_svc.fit(X_small_train, y_small_train)

print gs_svc.best_params_, gs_svc.best_score_ # {'C': 10.0, 'gamma': 0.001} 0.982

首先前面幾行就只是設定參數值而已, 通常這範圍跟經驗比較有關, 並沒有一定的測試數據

然後給定 training & test 的大小, 用切割好 train & test 設定 GridSearchCV

餵入 train & test, 他就會幫你找出給定的 data 中最適當的 C & gamma

通常我會測試 10-fold cross validation, 比較安心

那結果也找出來了, C = 10, gamma=1e-3, perfect, 上面也有執行過一次這樣的結果

相當方便吧~!!

Summary

就以一個剛拿到手的 dataset 想要用一個帶有變數的 learner 實驗怎麼開始比較好?

簡單來說就是可以先跑個 PCA 看看 class 的跟重要的 feature 之間分佈是否明顯

然後ShuffleSplit 切割好 data, 跟著跑 Cross validation 調參數, 找到之後再拿來使用

那個人覺得最好是上述整個流程都有一定程度熟悉, 再來直接使用 GridSearchCV

上面有交代不清或者是看不懂的相當歡迎提出來~(但是我有沒有辦法解答不一定XD)