scikit-learnでサポートベクトル回帰、及びそのパラメーター推計 with クロスバリデーションやってみる
サポートベクトル回帰(Support Vector Regression, SVR)の理論が大体脳内整理出来たので、実践もしたいぞと、そしてちょいとpythonを使う別件があるので、慣れの意味も込めてR言語ではなくpythonとその機械学習ライブラリであるscikit-learnを使ってやるぞとそういうことです。
scikit-learn自体のインストールはこの記事の最下部にある日本語のLINKを見れば良いと思う。
俺はpip使ってインストールしたような気がするけど、なにぶんずいぶんと昔なので忘れてしまった。pipで入れるなら
pip install scikit-learn
でOK。裏でコンパイルが走っていたような記憶があるので、C++のコンパイラいれておかないとだめかも。
windows用のバイナリファイルだと
にあるので、お手元のpythoのバージョンに合わせて入れちゃえばよさげ。
お試しの例題として用いるもの
ここでの例は
- -πからπまでの範囲のsin関数+noise
なデータ、つまりを満たすような(x,y)の組に対してサポートベクトル回帰をかますこととする。
データの作成法は、の区間から複数個一様にサンプリングした後、それをsin関数に食わせてノイズ(ここでは平均0,分散0.01ガウシアン)を乗せてyとした。
図示しておくと
なデータ(x,yのペア)が与えられているというシチュエーションを考えている。ここで言及した元データ＆グラフは以下のスクリプトで作成。
import numpy as np import matplotlib.pyplot as plt # X座標を適当にサンプリングして指定、それに合わせてy = sin(x) + noiseを生成 np.random.seed(1) x = np.sort(np.random.uniform(-np.pi, np.pi,100)) y = np.sin(x) + 0.1*np.random.normal(size=len(x)) # scikit-learnに突っ込むためにフォーマット変更 x = x.reshape((len(x),1)) # 全データの描画 plt.plot(x, y, 'o') plt.show()
サポートベクトル回帰(SVR)の実行
サポートベクトル回帰はscikit-learnにあるSVR関数を用いて実行する。上で作成した元々のデータを6:4の比率でトレーニング＆テストデータに分けて使用。
SVRの戻り値はsklearn.svm.classes.SVRクラスのオブジェクトになってて、こいつのpredictメソッドで予測、scoreメソッドで”予測のあてはまり度合い”を表す決定係数R^2を取得する事ができる*1
このオブジェクトについての詳細は
を見るべし。実際にトレーニングデータから学習した結果を使ってテストデータに対する当てはめを行ってみると
という感じになる(赤が予測で青が実際)。
コードは以下。
SVR関数にx,yを食わせる際、yをN次元のベクトルとした場合、xは1次元であろうとN×M次元の行列形式にしなければならない点に注意。
from sklearn import svm from sklearn import cross_validation #データを6割をトレーニング、4割をテスト用とする x_train, x_test, y_train, y_test = cross_validation.train_test_split(x, y, test_size=0.4) # 線でつないでplotする用にx_test・y_testをx_testの昇順に並び替える index = x_test.argsort(0).reshape(len(x_test)) x_test = x_test[index] y_test = y_test[index] # サポートベクトル回帰を学習データ使って作成 reg = svm.SVR(kernel='rbf', C=1).fit(x_train, y_train) # テストデータに対する予測結果のPLOT plt.plot(x_test, y_test, 'bo-', x_test, reg.predict(x_test), 'ro-') plt.show() # 決定係数R^2 print reg.score(x_test, y_test)
クロスバリデーション(Cross Validation, CV)
今手元にあるトレーニングデータを更にK個に分割し、そのうちの1個をテストデータ(A)、残りK-1個を改めてトレーニングデータ(B)とし、(B)を用いてモデルを構築し、(A)に対して予測を行い、モデルの妥当性を検証する事を考える。分割したデータがK個あるので、テストデータの選び方自体もK通りあり、そのK個に分割されたトレーニングデータそれぞれをテストデータとしてK回検証を行い、得られたK回の結果(スコア、何らかのモデルの良し悪しを決める指標）を平均して1つの評価指標の推定値を得ることが出来、これをK-foldクロスバリデーションといいます。その詳細は
を読めば十分わかるかなと。
scikit-learnでは"回帰"に関してモデルの良し悪しを決める指標として
- 決定係数R^2
- 平均二乗誤差(Mean Squared Error, MSE)
を指定出来る*2が、ここでMSEの値は”通常の(統計学|機械学習)”でいうMSEの値が-1倍されたものであり、これがスコアとして返ってくる点に気を付ける*3。
これをscikit-learn使って書くと以下のようになる。
# 5-foldクロスバリデーション。デフォルトはr^2がスコアとなって返ってくる。すなわち以下と等価 # cross_validation.cross_val_score(svm.SVR(), x, y, cv=5, scoring="r2") scores = cross_validation.cross_val_score(svm.SVR(), x, y, cv=5) # 5-foldクロスバリデーション毎のR^2の平均値とその±2σレンジ print("R^2(not adjusted): %0.2f (+/- %0.2f)" % (scores.mean(), scores.std() * 2)) #平均二乗誤差(MSE)をスコアに変更 scores = cross_validation.cross_val_score(svm.SVR(), x, y, cv=5, scoring="mean_squared_error") # 5-foldクロスバリデーション毎のMSEの平均値とその±2σレンジ print("MSE: %0.2f (+/- %0.2f)" % (-scores.mean(), scores.std() * 2))
クロスバリデーションとグリッドサーチを使ってSVRのパラメータを決める
単純にCVやるだけなら上でいいが、今まで所与だと思って放置していたモデルパラメータ、この例だと
- 罰則の強さを表す:C
- ガウシアンカーネルの広がりを表す:γ
を決めてやりたい。そのためにはGridSearchCV関数を用いる。
適当なパラメーターグリッド（組）に対してK-fold CV(ここでは5-fold）を繰り返し、スコア(ここではMSE)が最も良くなるパラメーターのSVRを出力する事ができる。
ちなみに以下のコードで
tuned_parameters = [{'kernel': ['rbf'], 'gamma': [10**i for i in range(-4,0)], 'C': [10**i for i in range(1,4)]}]
と書いている箇所を
tuned_parameters = [
{'kernel': ['rbf'], 'gamma': [10**i for i in range(-4,3)], 'C': [10**i for i in range(-3, 4)]},
{'kernel': ['linear'], 'C': [1, 10, 100, 1000]}
]
とすると、カーネルも同時に”パラメーター”として評価され、上記の通常言うようなパラメーターだけではなく、異なるカーネル間での結果比較も出来るのが素晴らしい*4。
このGridSearchCVを使って一番MSEがよくなるモデルを作成した結果が以下のグラフ。
「正答を青、最もMSEスコアの意味で(良|悪)いモデルの予測値を(赤|緑)」として示した。
コードは↓な感じ。一番スコアの良いモデルは
- best_estimator_
というプロパティに格納されているので簡単に取得できるが、（普通は用途のない）一番スコアの悪いモデルを出すのに苦労した。
from sklearn.grid_search import GridSearchCV #RBFカーネルのパラメーターγと罰則Cを複数個作ってその中で(スコアの意味で）良い物を探索(カーネルもパラメーターとして使用可能) tuned_parameters = [{'kernel': ['rbf'], 'gamma': [10**i for i in range(-4,0)], 'C': [10**i for i in range(1,4)]}] gscv = GridSearchCV(svm.SVR(), tuned_parameters, cv=5, scoring="mean_squared_error") gscv.fit(x_train, y_train) #一番スコア悪い&良い奴を出す params_min,_,_ = gscv.grid_scores_[np.argmin([x[1] for x in gscv.grid_scores_])] reg_min = svm.SVR(kernel=params_min['kernel'], C=params_min['C'], gamma=params_min['gamma']) reg_max = gscv.best_estimator_ #全トレーニングデータを使って再推計 reg_min.fit(x_train, y_train) reg_min.fit(x_train, y_train) #正答(青)＆良い(赤)＆悪い(緑)の結果をPLOT plt.plot(x_test, y_test, 'bo-',x_test, reg_max.predict(x_test), 'ro-',x_test, reg_min.predict(x_test), 'go-') plt.show()
その他メモ
- 大体の関数にn_jobsって引数があって、これを指定するだけで並列に走る数を制御できる。-1で全CPU使う
- グリッドサーチに関してはhttp://scikit-learn.org/stable/auto_examples/grid_search_digits.html#example-grid-search-digits-pyを見るのが良い
参考LINK
*1:なぜここでscoreメソッドが返す値として決定係数なのかはよくわかってない。SVRにおいて、決定係数があてはまりの良し悪しを判断するために使われる指標として普通なのだろうか???
*2:自作関数も可能。3.3. Model evaluation: quantifying the quality of predictions — scikit-learn 0.21.dev0 documentationを参照
*3:これはscikit-learnの仕様でmean_squared_error(MSE)関数をスコアとしてcross_val_score関数経由でCVやろうとすると、この関数内部で作られるscorerというスコア管理用のオブジェクトを作る際に使われるmake_scorer関数のgreater_is_better(値がでかい方がいい結果フラグ)引数が、MSEを使う場合false設定され、MSE自体の値が-1倍されるようになっている。make_scorer関数についてはsklearn.metrics.make_scorer — scikit-learn 0.19.2 documentationを参照
