Machine Learning with Scikit Learn (Part II)

Machine Learning with Scikit Learn (Part I)（2015/8/10）のつづき。今回は、後編のPartIIの動画の内容を簡単にまとめたい。

4.1 Cross Validation

ラベル付きデータが少ないときに有効な評価法であるK-fold cross-validationについての説明。訓練データをK個のサブセットに分割し、そのうち1つのサブセットをテストデータに残りK-1個のサブセットを訓練データにして評価する。これをテストデータを入れ替えながらK回評価し、その平均を求める。

この手順は下のように書ける。bool型のマスクを使ってテストデータと訓練データをわけている。

k = 5
n_samples = len(X)
fold_size = n_samples // k
scores = []
masks = []
for fold in range(k):
    test_mask = np.zeros(n_samples, dtype=bool)
    test_mask[fold * fold_size:(fold+1)*fold_size] = True
    X_test, y_test = X[test_mask], y[test_mask]
    X_train, y_train = X[~test_mask], y[~test_mask]
    classifier.fit(X_train, y_train)
    scores.append(classifier.score(X_test, y_test))
print np.mean(scores)

sklearn.cross_validationモジュールを使うと上と同じことが3行で書ける。

from sklearn.cross_validation import cross_val_score
scores = cross_val_score(classifier, X, y, cv=5)
print np.mean(scores)

データの分割方法は、StratifiedKFold、KFold、ShuffleSplit、LeaveOneOutなどが用意されている。たとえば、StratifiedKFoldのデータの分割方法は、

cv = StratifiedKFold(iris.target, n_folds=5)
for train, test in cv:
    print(test)

で確かめられる。可視化してみると

となる。行がCross Validationの各fold、列がirisの150個のデータである。青が訓練データ、赤がテストデータを意味している。irisデータは最初の50個がクラス1、次の50個がクラス2、次の50個がクラス3なので、各foldのテストデータはクラスの比率を維持したまま選ばれていることがわかる。

一方、単純なKFoldだと

と単純に等分にされる。データをシャッフルしないと各foldの正解ラベルが偏ってしまい、精度が大きくばらつくことになる。データをシャッフルするにはKFoldの引数にshuffle=Trueを入れておけばOK。

テストデータの選択をランダムで行うShuffleSplitもある。

ラベル付きデータ数が非常に少なく、K個のfoldに分割もできない場合は、テストデータとして1サンプルだけ選び、残りサンプルすべてを訓練データにするLeaveOneOutが使える。

cross_val_score()のデフォルトでは、StratifiedKFoldが使われるが、cvオプションに上記の任意のオブジェクトを指定できる。例えば、テストデータがランダムに選択されるShuffleSplitを使ってCross Validationをしたければ下のように書く。

cv = ShuffleSplit(len(iris.target), n_iter=5, test_size=.2)
cross_val_score(classifier, X, y, cv=cv)

4.2 Model Complexity and GridSearchCV

モデルのOverfitting/Underfittingとハイパーパラメータを最適化するグリッドサーチについての説明。通常の機械学習モデルは下のような関係が成り立つ。

引用:amueller/scipy_2015_sklearn_tutorial · GitHub

モデルが単純すぎるとUnderfittingぎみになり、訓練データとテストデータの両方とも精度が低い。一方、モデルを複雑化するとOverfittingぎみになり、訓練データの精度はすごく高くなるが、テストデータの精度が逆にすごく悪くなる。もっとも汎化された望ましいモデルは、テストデータの精度がもっとも高くなるSweet Spotの場所になる。

説明ではK近傍法の例が挙げられている。K近傍法のハイパーパラメータは近傍数のn_neighbors（K近傍法のKと同じ）である。このハイパーパラメータを変えて学習すると下のようになる。

一番左のK=2はモデルが複雑で訓練データの予測精度は高いがテストデータの精度は低いためOverfittingぎみ。一方、一番右のK=20はモデルが単純で訓練データに対してもテストデータに対しても精度は低いためUnderfittingぎみと解釈できる。真ん中のK=5が最適モデルとなる。

このような最適モデルを探す一般的な方法はなく、モデルのハイパーパラメータの組み合わせを力づくで評価して見つけるしかない。この最適なハイパーパラメータの探索にもCross Validationが使える。

from sklearn.cross_validation import cross_val_score, KFold
from sklearn.neighbors import KNeighborsRegressor

# generate toy dataset:
x = np.linspace(-3, 3, 100)
rng = np.random.RandomState(42)
y = np.sin(4 * x) + x + rng.normal(size=len(x))
X = x[:, np.newaxis]

cv = KFold(n=len(x), shuffle=True)

# for each parameter setting do cross_validation:
for n_neighbors in [1, 2, 3, 5, 10, 20]:
    scores = cross_val_score(KNeighborsRegressor(n_neighbors=n_neighbors), X, y, cv=cv)
    print "n_neighbors: %d, average score: %f" % (n_neighbors, np.mean(scores))

実行結果は、

n_neighbors: 1, average score: 0.586909
n_neighbors: 2, average score: 0.706250
n_neighbors: 3, average score: 0.738065
n_neighbors: 5, average score: 0.769667
n_neighbors: 10, average score: 0.740026
n_neighbors: 20, average score: 0.639899

となり、先のグラフからもわかるようにn_neighbors=5の平均精度がもっともよいことがわかる。validation_curve()を使うと横軸にハイパーパラメータ、縦軸に精度のグラフが簡単に描ける。

from sklearn.learning_curve import validation_curve
n_neighbors = [1, 2, 3, 5, 10, 20]
train_scores, test_scores = validation_curve(KNeighborsRegressor(), X, y, param_name="n_neighbors",
                                             param_range=n_neighbors, cv=cv)
plt.plot(n_neighbors, train_scores.mean(axis=1), label="train score")
plt.plot(n_neighbors, test_scores.mean(axis=1), label="test score")
plt.legend(loc="best")

このグラフからもn_neighbors=5でテストデータの精度が最大化されることがわかる。

K近傍法はハイパーパラメータが1つしかないが、複数ある場合も同じように探索できる。たとえば、SVMはCとgammaという2つのパラメータがあるため二重ループでパラメータのすべての組み合わせを評価する必要がある。

from sklearn.cross_validation import cross_val_score, KFold
from sklearn.svm import SVR

for C in [0.001, 0.01, 0.1, 1, 10]:
    for gamma in [0.001, 0.01, 0.1, 1]:
        scores = cross_val_score(SVR(C=C, gamma=gamma), X, y, cv=cv)
        print "C: %f, gamma: %f, average score: %f" % (C, gamma, np.mean(scores))

実行結果は、

C: 0.001000, gamma: 0.001000, average score: -0.006242
C: 0.001000, gamma: 0.010000, average score: -0.004179
C: 0.001000, gamma: 0.100000, average score: 0.004933
C: 0.001000, gamma: 1.000000, average score: 0.003641
C: 0.010000, gamma: 0.001000, average score: -0.003968
C: 0.010000, gamma: 0.010000, average score: 0.016451
C: 0.010000, gamma: 0.100000, average score: 0.098317
C: 0.010000, gamma: 1.000000, average score: 0.088333
C: 0.100000, gamma: 0.001000, average score: 0.018481
C: 0.100000, gamma: 0.010000, average score: 0.180360
C: 0.100000, gamma: 0.100000, average score: 0.526633
C: 0.100000, gamma: 1.000000, average score: 0.498146
C: 1.000000, gamma: 0.001000, average score: 0.194592
C: 1.000000, gamma: 0.010000, average score: 0.596778
C: 1.000000, gamma: 0.100000, average score: 0.672521
C: 1.000000, gamma: 1.000000, average score: 0.724113
C: 10.000000, gamma: 0.001000, average score: 0.593179
C: 10.000000, gamma: 0.010000, average score: 0.622648
C: 10.000000, gamma: 0.100000, average score: 0.669068
C: 10.000000, gamma: 1.000000, average score: 0.772627

となり、C=10とgamma=1.0のスコアが0.77と最大であることがわかる。

このようなハイパーパラメータの探索は、グリッドサーチと呼ばれる。scikit-learnではグリッドサーチを簡単に行うGridSearchCVが提供されている。

from sklearn.grid_search import GridSearchCV
param_grid = {'C': [0.001, 0.01, 0.1, 1, 10], 'gamma': [0.001, 0.01, 0.1, 1]}
grid = GridSearchCV(SVR(), param_grid=param_grid, cv=cv, verbose=3)

パラメータの組み合わせは辞書で与える。戻り値のgridに対してデータをfit()させるとハイパーパラメータの学習とモデル学習が順番に実行される。。

grid.fit(X, y)

最初に与えたデータを使ってCross Validationによる最適なハイパーパラメータ探索が行われ、その後に与えたデータすべてを使ってモデル学習が行われる。Cross Validationで見つかった最適なパラメータと最適な精度は下のコードで取得できる。

print grid.best_params_
print grid.best_score_

{'C': 10, 'gamma': 1}
0.772226773354

全体の手順をまとめると

1. 訓練データとテストデータに分割する
2. 訓練データを使ってCross Validationでモデルの最適なハイパーパラメータを見つける
3. 訓練データを使ってモデルを学習する
4. テストデータで精度を評価する

2と3はfit()でまとめて行われる。絵で描くと

引用:amueller/scipy_2015_sklearn_tutorial · GitHub

コードで書くと

from sklearn.cross_validation import train_test_split
X_train, X_test, y_train, y_test = train_test_split(X, y, random_state=0)

param_grid = {'C': [0.001, 0.01, 0.1, 1, 10],
              'gamma': [0.001, 0.01, 0.1, 1]}
cv = KFold(n=len(X_train), n_folds=10, shuffle=True)
grid = GridSearchCV(SVR(), param_grid=param_grid, cv=cv)
grid.fit(X_train, y_train)
grid.score(X_test, y_test)

4.3 Analyzing Model Capacity

使っているモデルがOverfittingなのかUnderfittingなのかを見極めるためには学習曲線（learning curve）を描くとよい。学習曲線は、横軸に訓練サンプル数、縦軸に精度（スコア）を描いたグラフである。scikit-learnでは、learning_curveメソッドで簡単に描ける。例では、SVMのカーネルにlinear、poly、rbfを使った場合の学習曲線を比較している。モデルの複雑さは、linear < poly < rbfとなる。

まずは、SVMのカーネルにlinearを使った場合。

from sklearn.learning_curve import learning_curve
from sklearn.svm import SVR
training_sizes, train_scores, test_scores = learning_curve(SVR(kernel='linear'),
                                                X, y, cv=10,
                                                scoring="mean_squared_error",
                                                train_sizes=[0.6, 0.7, 0.8, 0.9, 1.0])
plt.plot(training_sizes, train_scores.mean(axis=1), label="training scores")
plt.plot(training_sizes, test_scores.mean(axis=1), label="test scores")
plt.legend(loc="best")

学習曲線のスコアにmean squared errorを使うと戻り値はNegative MSEが返るので注意。MSEは小さいほうが良いがスコアは大きいほうがよいため符号を逆にして返す仕様らしい。learning_curve()の戻り値は訓練データ数のリスト、各訓練データ数で学習したときの訓練スコア、テストスコアとなる。10-fold Cross Validationで評価しているため訓練スコアとテストスコアの両方がまとめて求められる。

この学習曲線を見ると、訓練データ数を増やしても訓練スコア、テストスコアともに改善しない傾向があることがわかる。これはUnderfittingを意味する。つまり、そもそものモデルの表現力が貧弱なので頑張って訓練データ数を増やしても精度は改善しない。Underfittingなモデルに対しては、

• より複雑なモデルを使う（kernelをpolyやrbfにするなど）
• サンプルの特徴量を増やす
• 正則化を緩める（SVMのパラメータC）

などの対応を取るとよい。

次に、SVMのカーネルにrbfを使った場合。先のコードでkernel='rbf'に置き換えるだけ。

テストスコアに比べて訓練スコアが圧倒的に高い傾向がある。これはOverfittingを意味する。つまり、モデルが訓練データにフィットしすぎているため訓練データの精度は非常によくなるが、逆に汎化性能がなくなるためテストスコアが悪くなる。Overfittingなモデルに対しては、

• 特徴量を減らす
• より単純なモデルを使う
• 訓練データ数を増やす
• 正則化を強める

などの対応を取るとよい。

最後に、SVMのカーネルにpolyを使った場合。先のコードでkernel='poly'に置き換えるだけ。

訓練データ数を増やすと訓練スコアとテストスコアの精度の差が小さくなり、収束することがわかる（スコア軸のスケールが前の2つに比べて非常に小さいことに注目）。このような学習曲線が得られるとよいモデルであると判断できる。つまり、今回のデータではkernel=polyを使うのが最適ということがわかる。

4.4 Model Evaluation and Scoring Metrics

モデルのさまざまな評価法についての説明が続く。分類モデルのscore()はデフォルトでは精度（accuracy）、正解したラベルの割合を返す。

from sklearn.datasets import load_digits
from sklearn.cross_validation import train_test_split
from sklearn.svm import LinearSVC

digits = load_digits()
X, y = digits.data, digits.target
X_train, X_test, y_train, y_test = train_test_split(X, y)

classifier = LinearSVC().fit(X_train, y_train)
y_test_pred = classifier.predict(X_test)

print "Accuracy: %f" % classifier.score(X_test, y_test)

精度以外にもさまざまな評価法があり、sklearn.metricsに実装されている。多クラスの分類問題で有効な評価法が混同行列（confusion matrix）である。各クラスのサンプルがどんな間違われ方をしたのか直感的に把握できる。

from sklearn.metrics import confusion_matrix
confusion_matrix(y_test, y_test_pred)

[[53  0  0  0  0  0  0  0  0  0]
 [ 0 39  0  0  0  0  0  0  3  1]
 [ 0  1 39  0  0  0  0  0  0  0]
 [ 0  0  0 44  0  2  0  0  0  1]
 [ 0  0  0  0 40  0  0  0  0  1]
 [ 0  0  0  0  0 47  0  0  0  0]
 [ 0  0  0  0  1  0 45  0  0  0]
 [ 0  0  0  0  0  0  0 50  1  1]
 [ 0  3  0  0  0  0  0  0 39  2]
 [ 0  1  0  0  1  0  0  0  3 32]]

グラフで可視化するとより直感的にわかる。

plt.matshow(confusion_matrix(y_test, y_test_pred))
plt.colorbar()
plt.xlabel("Predicted label")
plt.ylabel("True label")

もう一つの便利な関数はclassification_report()。これは、各クラスの適合率（precision）、再現率（recall）、F1スコア（F1-score）、サポートを表形式で表示してくれる。

from sklearn.metrics import classification_report
print classification_report(y_test, y_test_pred)

             precision    recall  f1-score   support

          0       1.00      1.00      1.00        53
          1       0.89      0.91      0.90        43
          2       1.00      0.97      0.99        40
          3       1.00      0.94      0.97        47
          4       0.95      0.98      0.96        41
          5       0.96      1.00      0.98        47
          6       1.00      0.98      0.99        46
          7       1.00      0.96      0.98        52
          8       0.85      0.89      0.87        44
          9       0.84      0.86      0.85        37

avg / total       0.95      0.95      0.95       450

精度が役立たない状況として不均衡データ（imbalanced data）の例が挙げられている。不均衡データとはたとえば正例が10%で残りの90%が負例のように各クラスのサンプル数が非常に不均衡なデータ。こういうデータだとすべて負例と予測するだけで精度が90%になってしまう。

from sklearn.cross_validation import cross_val_score
from sklearn.svm import SVC

X, y = digits.data, digits.target == 3
print np.mean(cross_val_score(SVC(), X, y, cv=10))

この例では、数字の3が正例でそれ以外の数字が負例という二値分類問題をSVMで学習し、10-fold CVで評価している。予想通り精度は90%になる。すべてのテストデータに対して「3ではない」と推定するだけでよい。精度は役に立たない。昔これではまったことがあるのを思い出した（笑）

こういうときはROC曲線（ROC curve）とAUC（area under the curve）を使うとより正確に評価できる。ROC曲線は下のコードで描ける。

from sklearn.metrics import roc_curve, roc_auc_score

X, y = digits.data, digits.target == 3
X_train, X_test, y_train, y_test = train_test_split(X, y)

for gamma in [0.01, 0.05, 1]:
    plt.xlabel("False Positive Rate")
    plt.ylabel("True Positive Rate")
    svm = SVC(gamma=gamma).fit(X_train, y_train)
    decision_function = svm.decision_function(X_test)
    fpr, tpr, _ = roc_curve(y_test, decision_function)
    acc = svm.score(X_test, y_test)
    auc = roc_auc_score(y_test, svm.decision_function(X_test))
    plt.plot(fpr, tpr, label="acc:%0.2f auc:%0.2f" % (acc, auc), linewidth=3)
plt.legend(loc="best")

ROC曲線は横軸にFalse Positive Rate、縦軸にTrue Positive Rate（Recall）を割り当て、さまざまな閾値でこの2つ組をプロットしたグラフである。

横軸のFalse Positive Rate（FPR）は、本来は負例なのに誤って（False）正例（Positive）と判定してしまった割合なので小さいほどよい。一方、縦軸のTrue Positive Rate（TPR）は、正しく（True）正例（Positive）と判定できた割合なので大きい方がよい。つまり、カーブが左上に行くほど分類器の性能がよいことになる。このカーブの下の面積がAUCで1.0のとき最良となる。

この結果を見るとSVMのパラメータ（gamma）をいろいろ変えても先と同じように精度（accuracy）はすべて90%になってしまい優劣の比較ができない。しかし、ROC曲線とAUCの値はパラメータによって異なっており、gamma = 0.01（青の曲線）がもっともよいと判定できる。つまり、精度で評価するよりAUCで評価した方がよい場合がある。

ROC曲線は名前だけ聞いたことあったけどほとんど使ったことなかった。このページの説明はわかりやすいと思った。あとは朱鷺の杜のROC曲線のページも。

sklearnには上記以外にもさまざまなスコアがすでに実装されている。

from sklearn.metrics.scorer import SCORERS
print SCORERS.keys()

['f1', 'f1_weighted', 'f1_samples', 'recall_samples', 'recall_micro',
 'adjusted_rand_score', 'recall_macro', 'mean_absolute_error',
 'precision_macro', 'precision_weighted', 'f1_macro', 'recall_weighted',
 'accuracy', 'precision_samples', 'median_absolute_error', 'precision',
 'log_loss', 'precision_micro', 'average_precision', 'roc_auc', 'r2',
 'recall', 'mean_squared_error', 'f1_micro']

mean_squared_errorのようなよく使うものもあるけど、聞いたことないのもいろいろあるんだな。

ここまでで4章終わり。モデル評価は地味な感じがするけど、機械学習の応用では必須の知識になりそう。5章からはまた内容が一転してさまざまな分類器の詳細な説明が始まっている。長くなったのでつづきはPart IIIで。

Machine Learning with Scikit Learn (Part I)

今年の7月に開催されたSciPy2015の講演動画がEnthoughtのチャンネルで公開されている。今年も面白い講演が多いのでいろいろチェックしている。

今年の目標（2015/1/11）にPythonの機械学習ライブラリであるscikit-learnを使いこなすというのが入っているので、まずはscikit-learnのチュートリアルを一通り見ることにした。

Part IとPart IIを合わせると6時間以上あり非常に充実している。IPython Notebook形式の資料やデータは下記のGitHubアカウントで提供されている。ノートブックをダウンロードし、実際に手を動かしながらチュートリアルを進めると理解がより進むかもしれない。

あとで振り返りやすいように内容を簡単にまとめておきたい。

1.1 Introduction to Machine Learning

機械学習システムの流れ。教師あり学習と教師なし学習の違い。訓練データとテストデータは分けましょうという一般的なお話。

1.2 IPython Numpy and Matplotlib Refresher

scikit-learnのベースとなるNumPy、SciPy、Matplotlibのおさらい。言語処理では疎行列表現（Sparse Matrix）が効率的なのでscipy.sparseを使える。NumPy、SciPy、Matplotlibは非常に巨大なライブラリなので別途習得する必要がある。NumPyやMatplotlibについては以下の講演が過去にあった。これらもあとでチェックしておこう。

• Introduction to NumPy | SciPy 2015 Tutorial | Eric Jones - YouTube
• Matplotlib: past, present and future; SciPy 2013 Presentation - YouTube
• Anatomy of Matplotlib - Part 1 | SciPy 2014 | Benjamin Root - YouTube

np.newaxisを使った列ベクトルへの変換法は覚えておくと便利。1次元配列を行がサンプル、列が特徴量のscikit-learnのデータ形式に変換するときなどに使える。

# make into a column vector
print y[:, np.newaxis]

1.3 Data Representation for Machine Learning

scikit-learnのデータ表現は、行がサンプル、列が特徴量の2次元配列で表す。有名なデータセットはsklearn.datasetsモジュールでロードまたは生成できる。iris、digits、S-curve、olivetti_facesが例として取り上げられている。顔画像データなんてのも提供されているのか。あとで使ってみたいな。

1.4 Training and Testing Data

訓練データとテストデータの作り方。データはシャッフルされていないこともあるので必ずシャッフルする。データのシャッフルは、permutationをインデックスに指定すると簡単にできる。

import numpy as np
rng = np.random.RandomState(0)
permutation = rng.permutation(len(X))
X, y = X[permutation], y[permutation]

訓練データとテストデータの分割は、sklearn.cross_validation.train_test_split()を使うと簡単にできる。random_stateに乱数の種を指定すれば自分でシャッフルしなくてもOK。

from sklearn.cross_validation import train_test_split
train_X, test_X, train_y, test_y = train_test_split(X, y, train_size=0.5, random_state=1999)

2.1 Supervised Learning - Classification

教師あり学習の一種である分類について。sklearn.datasets.make_blobs()で人工的に生成した二次元データを使用。scikit-learnでは、どの機械学習アルゴリズムも

1. 分類器のオブジェクトを生成
2. fit()に訓練データを与えて分類器のパラメータを学習
3. predict()にテストデータを与えて予測
4. score()で分類器の性能を評価

というAPIで統一されているので覚えやすい。例としてロジスティック回帰とK-近傍法の例が紹介されている。

from sklearn.linear_model import LogisticRegression
classifier = LogisticRegression()
classifier.fit(X_train, y_train)
prediction = classifier.predict(X_test)
classifier.score(X_test, y_test)

学習したパラメータは、分類器オブジェクトのプロパティから取得できる。

print(classifier.coef_)
print(classifier.intercept_)

K-近傍法の例。ロジスティック回帰とメソッド名が同じなので覚えやすい。分類器を作成するときにアルゴリズムによって独自のパラメータを指定できる。

from sklearn.neighbors import KNeighborsClassifier
knn = KNeighborsClassifier(n_neighbors=1)
knn.fit(X_train, y_train)
knn.score(X_test, y_test)

2.2 Supervised Learning - Regression

教師あり学習の一種である回帰について。回帰は教師信号が連続値である点が分類と異なる点。例として線形回帰とK-近傍法が取り上げられている。K-近傍法は分類だけでなく、回帰にも使えるのか。

from sklearn.linear_model import LinearRegression
regressor = LinearRegression()
regressor.fit(X_train, y_train)

from sklearn.neighbors import KNeighborsRegressor
kneighbor_regression = KNeighborsRegressor(n_neighbors=1)
kneighbor_regression.fit(X_train, y_train)

2.3 Unsupervised Learning - Transformations and Dimensionality Reduction

教師なし学習のデータの標準化と次元削減について。

データの標準化を学習とは言わないような気がするけど、scikit-learnでは他の学習器と同じインタフェースで実装されている。たとえば、平均を0、分散を1にする正規化は、StandardScalarを使って以下のように書ける。

from sklearn.preprocessing import StandardScaler
scaler = StandardScaler()
scaler.fit(X)
X_scaled = scaler.transform(X)

標準化はデータの平均、標準偏差というパラメータを推定するという見方もできるためfit()というメソッドでも特に違和感はない。学習したモデルを使ってデータを変換する（標準化する）ときはpredict()ではなく、transform()というメソッドが定義されている。

主成分分析もほとんど同じように書ける。

from sklearn.decomposition import PCA
pca = PCA(n_components=2)
pca.fit(X_blob)
X_pca = pca.transform(X_blob)

n_componentsで固有値が大きい順にいくつの固有ベクトル（軸）を取るか指定している。この例では、n_components=2なのでデータを2次元平面に圧縮することを意味する。次元削減はtransform()で行える。X_blobが元のデータでX_pcaが低次元空間に写像したデータになる。8x8ピクセル=64次元の数字画像データを2次元に写像すると下の結果が得られる。

高次元データを低次元のデータに縮約する手法をまとめて多様体学習（manifold learning）と呼ぶようだ。有名どころは主成分分析や多次元尺度構成法（MDS）だが、それ以外にも様々な手法が提案されており、一部はsklearn.manifoldに実装されている。たとえば、t-SNEを同じ数字画像データに適用すると

となり2次元空間で元の数字をきれいに分類できるようになる。この分野は非常に興味があるので勉強してみたいところ。

2.4 Unsupervised Learning - Clustering

教師なし学習の代表であるクラスタリングについて。有名なK-meansが例として取り上げられている。

from sklearn.cluster import KMeans
kmeans = KMeans(n_clusters=3, random_state=42)
labels = kmeans.fit_predict(X)

K-meansはユーザがクラスタ数をあらかじめ指定する必要がある。この例ではクラスタ数を3としている。講演では、K-means以外の様々なクラスタリング手法も紹介されている。

こんなにいろいろあったのか。K-meansとWard法くらいしか知らなかった。私の主観的なクラスタの判断だとDBSCANってのが一番近いかな？ここもしっかり勉強してみたいところ。個人的に教師あり学習よりも教師なし学習や表現学習の分野の方が興味ある。

2.5 Review of Scikit-learn API

scikit-learnで使われる fit(), predict(), transform() などのメソッド名についてのまとめ。教師あり学習、教師なし学習で提供されているメソッド名が異なる。さらにアルゴリズムによってもメソッドの意味が異なる場合がある。たとえば、score()は教師あり学習では分類精度、平均二乗誤差を返すが、一部の教師なし学習では尤度を返す。このノートブックにまとめられているので適宜参照したい。

3.1 Case Study - Supervised Classification of Handwritten Digits

ここからしばらく具体的なケーススタディが続く。最初は手書き数字認識。

1. 手書き数字データはsklearn.datasets.load_digits()でロードできる。
2. データ数が多いときはPCAより高速なsklearn.decomposition.RandomizedPCAを使ってデータの分布を可視化する。
3. PCAは線形の次元圧縮なのでデータの面白い関係を見逃している可能性がある。そのような場合は、sklearn.manifold.Isomapのような非線形な多様体学習手法を使う。
4. まずは簡単で高速な分類器（たとえば、sklearn.naive_bayes.GaussianNB）を使ってベースラインとする。
5. テストデータを用いて定量的な評価をする。分類問題ではmodel.score()が使える。また、適合率・再現率・F1値をまとめて返すsklearn.metrics.classification_report()やConfusion Matrixを返すsklearn.metrics.confusion_matrix()も便利。

from sklearn import metrics
print metrics.classification_report(expected, predicted)

             precision    recall  f1-score   support

          0       1.00      1.00      1.00        37
          1       0.68      0.91      0.78        43
          2       1.00      0.39      0.56        44
          3       0.76      0.84      0.80        45
          4       0.93      0.71      0.81        38
          5       0.98      0.90      0.93        48
          6       0.95      1.00      0.97        52
          7       0.70      0.98      0.82        48
          8       0.55      0.77      0.64        48
          9       1.00      0.57      0.73        47

avg / total       0.85      0.81      0.80       450

print metrics.confusion_matrix(expected, predicted)

[[37  0  0  0  0  0  0  0  0  0]
 [ 0 39  0  0  0  0  1  0  3  0]
 [ 0  9 17  3  0  0  0  0 15  0]
 [ 0  0  0 38  0  0  0  2  5  0]
 [ 0  1  0  0 27  0  2  8  0  0]
 [ 0  1  0  1  0 43  0  3  0  0]
 [ 0  0  0  0  0  0 52  0  0  0]
 [ 0  0  0  0  1  0  0 47  0  0]
 [ 0  5  0  1  0  1  0  4 37  0]
 [ 0  2  0  7  1  0  0  3  7 27]]

3.2 Methods - Unsupervised Preprocessing

Labeled Faces in the Wildの顔画像データに対して主成分分析を適用する。いわゆる固有顔（eigenface）の実験。このデータもメソッド一つでロードできる。

from sklearn import datasets
lfw_people = datasets.fetch_lfw_people(min_faces_per_person=70, resize=0.4, data_home='datasets')

これが訓練データの平均顔。著作権の問題からブッシュ元大統領の顔写真が多いんだってさ（笑）

これが各固有ベクトルを画像化した固有顔。

最初の数個の軸は画像の明るさの違いを表していることがわかる。これら固有顔の線形結合でいろいろな顔が表せる。

3.3 Case Study - Face Recognition with Eigenfaces

オリジナル画像の1850次元の顔画像データを150次元の固有顔空間に次元圧縮したデータを用いて顔認識の実験。分類にはサポートベクトルマシン（SVM）を使用。

from sklearn import svm
clf = svm.SVC(C=5., gamma=0.001)
clf.fit(X_train_pca, y_train)

メソッド名はsklearn.svm.SVCという名前なので注意。Support Vector Classificationの略。sklearn.svm.SVRという回帰版もある。こちらはSupport Vector Regressionの略。

オリジナル画像データをそのまま使うと40%くらいの精度しか出ないが、固有顔に次元圧縮してから分類すると85%くらいの精度が出る。

PCAをしてからSVMのような流れ作業はPipelineの機能を使うと便利。

from sklearn.pipeline import Pipeline
from sklearn import metrics
clf = Pipeline([('pca', decomposition.RandomizedPCA(n_components=150, whiten=True, random_state=1999)),
                ('svm', svm.SVC(C=5., gamma=0.001))])
clf.fit(X_train, y_train)
y_pred = clf.predict(X_test)
print metrics.classification_report(y_pred, y_test)

3.4 Methods - Text Feature Extraction

ここから対象がテキストになる。英語のテキストを特徴量に変換するbag-of-wordsの作り方について。CountVectorizerを使うとテキスト集合を単語の出現頻度を用いたbag-of-words表現に簡単に変換できる。

from sklearn.feature_extraction.text import CountVectorizer
X = ["Some say the world will end in fire,",
     "Some say in ice."]
vectorizer = CountVectorizer()
vectorizer.fit(X)
X_bag_of_words = vectorizer.transform(X).toarray()

単語の頻度TFだけでなく、逆文書頻度IDFも考慮したTF-IDFもTfidfVectorizerを使うと簡単に求められる。IDFを考慮すると多くの文書に出てくる単語が一般的な単語と判断されてその重要度が低くなる。

from sklearn.feature_extraction.text import TfidfVectorizer
tfidf_vectorizer = TfidfVectorizer()
tfidf_vectorizer.fit(X)
tfidf = tfidf_vectorizer.transform(X).to_array()

単語単位ではなく、N-gram単位でbag-of-wordsやTF-IDFを求めたい場合は、ngram_rangeオプションを指定すればOK。

# look at sequences of tokens of minimum length 2 and maximum length 2
bigram_vectorizer = CountVectorizer(ngram_range=(2, 2))
bigram_vectorizer.fit(X)

日本語では形態素解析が必要なのでこれらのメソッドはそのままでは使えなさそう。誰か日本語用の拡張も書いてたりするかな？

3.5 Case Study - SMS Spam Detection

テキストをbag-of-words表現に変換し、スパムのテキスト分類を行う例が紹介されている。訓練データの文章数は4180でボキャブラリ数は7464。ロジスティック回帰を使うと98%の分類精度が出ている。

3.6 Case Study - Titanic Survival

講演では時間が足りなくて取り上げられなかった（PartIIの最後に時間があまったため説明されていた）。カテゴリ特徴量を使うときは、各カテゴリに数値を割り当てるのではなく、いわゆる1-of-K表現にしましょうというお話。DictCategorizerを使うと簡単にできる。以下の例では、カテゴリ特徴量 city = {Dubai, London, San Francisco}と数値特徴量 temperature が混ざったデータをベクトル化している。

measurements = [
    {'city': 'Dubai', 'temperature': 33.},
    {'city': 'London', 'temperature': 12.},
    {'city': 'San Francisco', 'temperature': 18.},
]
from sklearn.feature_extraction import DictVectorizer
vec = DictVectorizer()
vec.fit_transform(measurements).toarray()

実行結果は、

array([[  1.,   0.,   0.,  33.],
       [  0.,   1.,   0.,  12.],
       [  0.,   0.,   1.,  18.]])

カテゴリ特徴量のみ1-of-K表現になって、数値特徴量はそのままになることがわかる。より複雑なデータセットとしてカテゴリ特徴量と数値特徴量が混ざったデータであるTitanic Dataを例に説明している。

ここで動画のPart Iは終わり。Part IIにつづきます。

Theanoによる畳み込みニューラルネットワークの実装 (2)

Theanoによる畳み込みニューラルネットワークの実装 (1)（2015/6/26）のつづき。今回は前回できなかった

• ConvLayerとPoolingLayerの分離
• ReLUの導入

を試してみた。

ConvLayerとPoolingLayerの分離

Deep Learning Tutorialの実装では、LeNetConvPoolLayer()というクラスで畳み込みとプーリングがセットで行われていた。

class LeNetConvPoolLayer(object):
    """畳み込みニューラルネットの畳み込み層＋プーリング層"""
    def __init__(self, rng, input, image_shape, filter_shape, poolsize=(2, 2)):

        # 入力の特徴マップとフィルタの畳み込み
        conv_out = conv.conv2d(
            input=input,
            filters=self.W,
            filter_shape=filter_shape,
            image_shape=image_shape)

        # Max-poolingを用いて各特徴マップをダウンサンプリング
        pooled_out = downsample.max_pool_2d(
            input=conv_out,
            ds=poolsize,
            ignore_border=True)

        # バイアスを加える
        self.output = T.tanh(pooled_out + self.b.dimshuffle('x', 0, 'x', 'x'))

いろいろ論文を読むと畳み込みとプーリングは必ずしもセットで行う必要はなく、畳み込みを数回繰り返した後にプーリングを1回という実装もありえるみたい。このクラスのままではそのようなネットワーク構成が作れないのでまずは畳み込み層とプーリング層を分離することにした。

また、上の実装では、

1. 畳み込み層
2. プーリング層
3. バイアスを加える
4. 活性化関数（tanh）

という順番になっていたが、人工知能学会の特集論文では、

1. 畳み込み層
2. バイアスを加える
3. 活性化関数
4. プーリング層

と書いてあった。この場合、学習パラメータである重みWとバイアスbは全部畳み込み層のクラスに収まり、プーリング層には学習するパラメータがなくなる。畳み込み層とプーリング層を分離する実装ではこっちの方が都合がよさそうなので採用した。

まずは畳み込み層のクラス。

class ConvLayer(object):
    """畳み込みニューラルネットの畳み込み層"""
    def __init__(self, rng, input, image_shape, filter_shape):
        # 入力の特徴マップ数は一致する必要がある
        assert image_shape[1] == filter_shape[1]

        fan_in = np.prod(filter_shape[1:])
        fan_out = filter_shape[0] * np.prod(filter_shape[2:])

        W_bound = np.sqrt(6.0 / (fan_in + fan_out))
        self.W = theano.shared(
            np.asarray(rng.uniform(low=-W_bound, high=W_bound, size=filter_shape),
                       dtype=theano.config.floatX),
            borrow=True)

        b_values = np.zeros((filter_shape[0],), dtype=theano.config.floatX)
        self.b = theano.shared(value=b_values, borrow=T)

        # 入力の特徴マップとフィルタの畳み込み
        conv_out = conv.conv2d(
            input=input,
            filters=self.W,
            filter_shape=filter_shape,
            image_shape=image_shape)

        # バイアスを加える
        self.output = T.tanh(conv_out + self.b.dimshuffle('x', 0, 'x', 'x'))

        self.params = [self.W, self.b]

畳み込み層にバイアスを加える処理と活性化関数（tanh）を通す処理を書いている。パラメータはWとbの両方が含まれる。次にプーリング層のクラス。

class PoolingLayer(object):
    """畳み込みニューラルネットのプーリング層
    この実装ではプーリング層にパラメータはない"""
    def __init__(self, rng, input, poolsize=(2, 2)):
        # Max-poolingを用いて各特徴マップをダウンサンプリング
        pooled_out = downsample.max_pool_2d(
            input=input,
            ds=poolsize,
            ignore_border=True)

        self.output = pooled_out

こちらには学習すべきパラメータはないので、self.paramsは定義しなくていい。

畳み込み層とプーリング層を分離したので見た目のレイヤの数は多くなる。前回と同じ構成の畳み込みニューラルネットを構築するには下のように書く必要がある。ConvLayerとPoolingLayerを分離したので少し大変。

    # 入力
    # 入力のサイズを4Dテンソルに変換
    # batch_sizeは訓練画像の枚数
    # チャンネル数は1
    # (28, 28)はMNISTの画像サイズ
    layer0_input = x.reshape((batch_size, 1, 28, 28))

    layer0 = ConvLayer(rng,
                input=layer0_input,
                image_shape=(batch_size, 1, 28, 28),
                filter_shape=(20, 1, 5, 5))

    layer1 = PoolingLayer(rng,
                          input=layer0.output,
                          poolsize=(2, 2))

    layer2 = ConvLayer(rng,
                       input=layer1.output,
                       image_shape=(batch_size, 20, 12, 12),
                       filter_shape=(50, 20, 5, 5))

    layer3 = PoolingLayer(rng,
                          input=layer2.output,
                          poolsize=(2, 2))

    # 隠れ層への入力
    layer4_input = layer3.output.flatten(2)

    # 全結合された隠れ層
    layer4 = HiddenLayer(rng,
        input=layer4_input,
        n_in=50 * 4 * 4,
        n_out=500,
        activation=T.tanh)

    # 最終的な数字分類を行うsoftmax層
    layer5 = LogisticRegression(input=layer4.output, n_in=500, n_out=10)

また、パラメータは畳み込み層（layer0とlayer2）、隠れ層（layer4）、softmax層（layer5）にしかないので

    # パラメータ
    params = layer5.params + layer4.params + layer2.params + layer0.params

となる。GPUを使って学習してみると

Using gpu device 0: GeForce GTX 760 Ti OEM
Optimization complete.
Best validation score of 0.950000 % obtained at iteration 11700, with test performance 1.060000 %
The code for file convnet2.py ran for 42.09m

となった。学習時間は41分から43分と少し遅くなった。テストエラー率は0.93%から1.06%になったのでこちらも少し悪化している。2回繰り返してみたけどどちらも同じ傾向だったのでもしかしたらDeep Learning Tutorialの順番の方がよかったのかもしれない。

ReLUの導入

ニューラルネットの非線形な活性化関数には、sigmoidやtanhが使われることが多かったが、近年ではRectified Linear Unit、略してReLUを使うことが多いそうだ。

Rectifier (neural networks) - Wikipedia, the free encyclopedia

下のように定義が非常に単純。ReLUを使うと勾配が減衰しにくくなるため収束性や学習速度が向上するという。

というわけで活性化関数をtanhからReLUに変えて学習させてみよう。ReLUの実装は非常に簡単。Theanoの機能を使うと

def relu(x):
    """Rectified Linear Unit"""
    return T.switch(x < 0, 0, x)

と書ける。Tはtheano.tensorの省略形。xはTensorVariableなので普通のif文では書けなかった。次に活性化関数にtanh()を使っていたところを全部reluに置換する。

    class ConvLayer(object):
        """畳み込みニューラルネットの畳み込み層"""
        def __init__(self, rng, input, image_shape, filter_shape):
            ...

            # バイアスを加える
            self.output = relu(conv_out + self.b.dimshuffle('x', 0, 'x', 'x'))

    ....

    # 全結合された隠れ層
    layer4 = HiddenLayer(rng,
        input=layer4_input,
        n_in=50 * 4 * 4,
        n_out=500,
        activation=relu)

同じくGPUで学習してみると

Using gpu device 0: GeForce GTX 760 Ti OEM
Optimization complete.
Best validation score of 1.200000 % obtained at iteration 15900, with test performance 1.060000 %
The code for file convnet2.py ran for 42.57m

となった。えー、最終的な学習時間や精度は全然改善していないようにみえる。もしかして収束が速いのかな？sigmoidとreluのエラー率の推移をグラフ化してみよう。

まあ、多少ReLUの方がtanhより収束が速いけど巷で噂になっているほどほどすごいわけでもないような・・・ネットワーク構成を変えるともっと効いたりするのかな？