今さら新しみのないサポートベクターマシン(SVM)のサンプルコード（Python,Scikit-learn）

・2年くらい前は、ネット上にあんまり日本語ドキュメントがなかったですが、今はたくさんあるサポートベクターマシン(Support Vector Machine,SVM)の回帰分析ver(SVR)をパッケージから呼び出すコードと実行結果を載せます。

f:id:russENG:20160320090508p:plain

１．SVMとは？
２．例題
３．いかに正しい予測器を選べば良いか
４．コードと結果

１．SVMとは？

　ちゃんとした説明は他のドキュメントで読んでもらえばよいとして、ざっくりと。
　
　※理解しなくても使えます（し、私も理解していない。）
　
　・１個の「パーセプトロン」を最適化問題として定式化して解く。この際の「マージン最大化」というアイディアがSVMの独自性。
　・多くの場合、単純なベクトル空間ではなく「カーネル誘導特徴空間」で判別・回帰作業を行なう。これにより、モデルの表現力が増し、SVMは爆発的な流行を見せた（ディープラーニングが流行る前）。
　・「パーセプトロン」とはもっとも素朴な識別機（私もこの辺りまでなら数学的に追える）であり、これをたくさん使うとk-NN法、直列に重ね合わせて、誤差逆伝播法というアイディアを実装したのが、いわゆるニューラルネットワーク。

　・最適化問題の適用時のいろんな工夫により、次元の呪いにかからない。また、局所最適解に陥らない。そしてスパース行列（疎行列）を用いるので計算が速い。

以下、私感やメモ。

　　・かなり表現力のあるモデルなので、過学習に注意する機構（多くの場合、交差検定（cross-validation）によるパラメータ制御）を使わないと全く使えない。
　　・訓練データ数が少なすぎると使えない。
　　・説明変数－被説明変数に相関関係がある場合はもちろん、非線形な関係であったとしても同じ枠組みで簡単に扱えるので便利。とりあえずSVMにぶち込んで結果がよければ、「何かしらの関係がある」ということになるので、有用なデータなのかどうかの判別機にもなるのかも。
　　・もちろん機械学習ツールなので、変数間の関係を学習しているだけ。変数の確率分布や、どの変数が効いているかとかを可視化することはできない。説明変数が多くてただぶち込んでいるだけの場合、あまり現象の理解には繋がっていない。
　　・ネット上を見ていると、こういうパッケージを使うのでなくて、自分で式から実装している例もずいぶんと増えて、みんな頭良いなあ、と思う。実装する場合、ラグランジュの未定乗数法という最適化数学の勉強にもなるのでカッコイイ。
　

２．例題

　ボストンの住宅価格データ(Boston House Prices Dataset)というデータが、scikit-learnのデータセットの中に。ざっくりと英語を呼んだところ、街ごとの統計データからその町の住宅価格（中央値）を当てる、というもの。部屋数の平均や、河川が近いかどうか、黒人が多いかどうか、築年数の平均なんかが説明変数として入っている。詳しくは後で記事にするかも。

３．いかに正しい予測器を選べばよいか

scikit-learnのドキュメント内でこんなの見つけた(今更)。参考になれば。
f:id:russENG:20160311074949p:plain *1

これに従うと、回帰問題の中で「説明変数のうち、被説明変数と関係あるものが少ない」とあらかじめわかっている場合は、Lassoなりを使ってくださいよ、と。で、LinearSVRから試せよ*2、ということです。

今回は、そんなの知らんがな(´・ω・｀)　でとりあえずSVR(しかも線形じゃないやつ)を使います。

タイトルにもあるけれど、機械学習の中での話なので、一般的な統計はまた別のお話（として私は頭の中で整理している。）

４．コードと結果

# -*- coding: utf-8 -*-
#import numpy
import matplotlib.pyplot as plt
from sklearn import datasets
from sklearn import svm
from sklearn import cross_validation
from sklearn.grid_search import GridSearchCV
from sklearn import preprocessing


# Bostonデータの読み込み
boston = datasets.load_boston()
X = boston.data
Y = boston.target

# データの6割をキャリブレーション用(校正曲線作成用)、4割をバリデーション用（校正曲線検証用）とする
# ※cross_validationというモジュールを使っているが、データのランダムな分割に使っているだけ。
X_cali, X_vali, Y_cali, Y_vali = cross_validation.train_test_split(X, Y, test_size=0.4)

# 説明変数の標準化
scaler = preprocessing.StandardScaler().fit(X_cali)
X_cali = scaler.transform(X_cali)
X_vali = scaler.transform(X_vali)

# 交差検定により、適切なパラメータを探索させる（グリッドサーチ）
# デフォルト：γ∈2^-15~2^3 (-15,4) C∈2^-5~2^15(-5,16) (出典有り)


# カーネルにRBFを適用
tuned_parameters = [{'kernel': ['rbf'], 'gamma': [2**i for i in range(-15,4)], 'C': [2**i for i in range(-5,16)]}]
gscv = GridSearchCV(svm.SVR(), tuned_parameters, cv=10,)
gscv.fit(X_cali, Y_cali)

reg_max = gscv.best_estimator_ #交差検定の結果、もっとも成績のよい判別器

# キャリブレーション結果の出力
plt.subplot(1,2,1)
plt.scatter(reg_max.predict(X_cali), Y_cali)
plt.title("calibration")
ymax = max(Y)
ymin = min(Y)
plt.xlim(ymin, ymax)
plt.ylim(ymin, ymax)
plt.xlabel('Actual') # x軸のラベル
plt.ylabel('Predicted') # y軸のラベル
scocali = "score = " + str(reg_max.score(X_cali, Y_cali))
plt.text(ymin,ymax*0.8, scocali)#位置は適当なので状況に応じて直すべき

# バリデーション結果の出力
plt.subplot(1,2,2)
plt.scatter(reg_max.predict(X_vali), Y_vali)
plt.title("validation")
plt.xlim(ymin, ymax)
plt.ylim(ymin, ymax)
plt.xlabel('Actual') 
plt.ylabel('Predicted') 
scovali = "score = " + str(reg_max.score(X_vali, Y_vali))
plt.text(ymin,ymax*0.8, scovali)
plt.show()