無料で使えるシステムトレードフレームワーク「Jiji」 をリリースしました!

・OANDA Trade APIを利用した、オープンソースのシステムトレードフレームワークです。
・自分だけの取引アルゴリズムで、誰でも、いますぐ、かんたんに、自動取引を開始できます。

機械学習手習い: サポートベクターマシン

「入門 機械学習」手習い、12日目。「12章 モデル比較」です。

www.amazon.co.jp

最後のアルゴリズムサポートベクターマシン(SVM)を学び、最後に同じデータセットにロジスティック回帰やk近傍法など、今まで学んできたアルゴリズムを適用して比較します。

# 前準備
> setwd("12-Model_Comparison/")

サポートベクターマシン(SVM)

サポートベクターマシンは、分類モデルの一つで、ロジスティック回帰と違い、非線形の決定境界を持つデータもうまく分類することができます。

例えば、以下のようなデータ。決定境界が1つの線で表せないため、ロジスティック回帰ではうまく分類できません。

> library('ggplot2')

# データを読み込み
> df <- read.csv(file.path('data', 'df.csv'))
> head(df)
          X          Y Label
1 0.2655087 0.52601906     1
2 0.3721239 0.07333542     1
3 0.5728534 0.84974175     1
4 0.9082078 0.42305801     0

> ggplot(df, aes(x = X, y = Y, color = factor(Label))) + geom_point()
> ggsave(filename="plot01.png")

f:id:unageanu:20160121150757p:plain

試しに、ロジスティック回帰を使って、任意の点がどのラベルに属するか判定してみます。

> logit.fit <- glm(Label ~ X + Y,
    family = binomial(link = 'logit'), data = df)

> logit.predictions <- ifelse(predict(logit.fit) > 0, 1, 0)
> mean(with(df, logit.predictions == Label)) 
[1] 0.5156

精度は51.5%。クラスは2つかないので、ランダムに選んだ場合とほぼ変わらない精度です。

次に、SVMで分類した場合の精度も計測してみます。

> library('e1071')

> svm.fit <- svm(Label ~ X + Y, data = df)
> svm.predictions <- ifelse(predict(svm.fit) > 0, 1, 0)
> mean(with(df, svm.predictions == Label))
[1] 0.7204

72%になりました。ロジスティック回帰よりうまく分類できているようです。

各モデルがどのような判定を行っているか、図示してみます。

> library("reshape")
> df <- cbind(df, data.frame(Logit = ifelse(predict(logit.fit) > 0, 1, 0),
  SVM = ifelse(predict(svm.fit) > 0, 1, 0)))

> predictions <- melt(df, id.vars = c('X', 'Y'))

> ggplot(predictions, aes(x = X, y = Y, color = factor(value))) +
  geom_point() + facet_grid(variable ~ .)
> ggsave(filename="plot02.png")

f:id:unageanu:20160121150758p:plain

上から、分類対象のデータ、ロジスティック回帰で各点を分類した結果、SVMで各点を分類した結果、です。

ロジスティック回帰は全部0にしている・・・。一方、SVMだと元データの特徴を完全ではないですがうまく捉えているようです。

カーネルトリック

SVMでは、カーネルトリックという手法を使って、非線形の決定境界を生成できるようになっています。 カーネルトリックは、数学的な変換で元のデータセットを新しい空間に移動することで、決定境界を線形でも記述しやすくすします。

何種類かのカーネルで分類を試し、結果を比較してみます。 svm 関数では、引数で利用するカーネルを変更できるので、それを使います。

> df <- df[, c('X', 'Y', 'Label')]

# 線形カーネル
> linear.svm.fit <- svm(Label ~ X + Y, data = df, kernel = 'linear')

# 多項式カーネル
> polynomial.svm.fit <- svm(Label ~ X + Y, data = df, kernel = 'polynomial')

# ガウスカーネル
> radial.svm.fit <- svm(Label ~ X + Y, data = df, kernel = 'radial')

# シグモイドカーネル
> sigmoid.svm.fit <- svm(Label ~ X + Y, data = df, kernel = 'sigmoid')

> df <- cbind(df,
            data.frame(LinearSVM = ifelse(predict(linear.svm.fit) > 0, 1, 0),
                       PolynomialSVM = ifelse(predict(polynomial.svm.fit) > 0, 1, 0),
                       RadialSVM = ifelse(predict(radial.svm.fit) > 0, 1, 0),
                       SigmoidSVM = ifelse(predict(sigmoid.svm.fit) > 0, 1, 0)))

> predictions <- melt(df, id.vars = c('X', 'Y'))

> ggplot(predictions, aes(x = X, y = Y, color = factor(value))) +
  geom_point() + facet_grid(variable ~ .)
> ggsave(filename="plot03.png")

f:id:unageanu:20160121150759p:plain

線形カーネル多項式カーネルはロジスティック回帰と同様、うまく分類できていない感じ。 ガウスカーネルは、正解に近い境界を生成できています。

今度は、パラメータを変えたパターンを試してみます。 多項式カーネルでは、次数をパラメータとして指定できるので、3,5,10,12のパターンで判定を行ってみます。

# 次数を変更したパターン
> polynomial.degree3.svm.fit <- svm(
  Label ~ X + Y, data = df, kernel = 'polynomial', degree = 3)
> polynomial.degree5.svm.fit <- svm(
  Label ~ X + Y, data = df, kernel = 'polynomial', degree = 5)
> polynomial.degree10.svm.fit <- svm(
  Label ~ X + Y, data = df, kernel = 'polynomial', degree = 10)
> polynomial.degree12.svm.fit <- svm(
  Label ~ X + Y, data = df, kernel = 'polynomial', degree = 12)

> df <- df[, c('X', 'Y', 'Label')]
> df <- cbind(df, data.frame(
  Degree3SVM = ifelse(predict(polynomial.degree3.svm.fit) > 0, 1, 0),
  Degree5SVM = ifelse(predict(polynomial.degree5.svm.fit) > 0, 1, 0),
  Degree10SVM = ifelse(predict(polynomial.degree10.svm.fit) > 0, 1, 0),
  Degree12SVM = ifelse(predict(polynomial.degree12.svm.fit) > 0, 1, 0)
))
> predictions <- melt(df, id.vars = c('X', 'Y'))
> ggplot(predictions, aes(x = X, y = Y, color = factor(value))) +
  geom_point() + facet_grid(variable ~ .)
> ggsave(filename="plot04.png")

f:id:unageanu:20160121150800p:plain

次数が3,5の場合は、指定なしの場合と変わりませんが、10,12だと少し判定できるようになっています。

次は、ガウスカーネルcost パラメータを変えてみます。 cost正則化の強さを示すパラメータで、値を大きくすると正則化が強く働き、訓練データに当てはまりにくくなります。

> radial.cost1.svm.fit <- svm(Label ~ X + Y, data = df, kernel = 'radial', cost = 1)
> radial.cost2.svm.fit <- svm(Label ~ X + Y, data = df, kernel = 'radial', cost = 2)
> radial.cost3.svm.fit <- svm(Label ~ X + Y, data = df, kernel = 'radial', cost = 3)
> radial.cost4.svm.fit <- svm(Label ~ X + Y, data = df, kernel = 'radial', cost = 4)

> df <- df[, c('X', 'Y', 'Label')]
> df <- cbind(df,
            data.frame(Cost1SVM = ifelse(predict(radial.cost1.svm.fit) > 0, 1, 0),
                       Cost2SVM = ifelse(predict(radial.cost2.svm.fit) > 0, 1, 0),
                       Cost3SVM = ifelse(predict(radial.cost3.svm.fit) > 0, 1, 0),
                       Cost4SVM = ifelse(predict(radial.cost4.svm.fit) > 0, 1, 0)))

> predictions <- melt(df, id.vars = c('X', 'Y'))
> ggplot(predictions, aes(x = X, y = Y, color = factor(value))) +
  geom_point() + facet_grid(variable ~ .)
> ggsave(filename="plot05.png")

f:id:unageanu:20160121150801p:plain

最後に、 シグモイドカーネルgamma パラメータを試して終わり。

> sigmoid.gamma1.svm.fit <- svm(Label ~ X + Y, data = df, kernel = 'sigmoid', gamma = 1)
> sigmoid.gamma2.svm.fit <- svm(Label ~ X + Y, data = df, kernel = 'sigmoid', gamma = 2)
> sigmoid.gamma3.svm.fit <- svm(Label ~ X + Y, data = df, kernel = 'sigmoid', gamma = 3)
> sigmoid.gamma4.svm.fit <- svm(Label ~ X + Y, data = df, kernel = 'sigmoid', gamma = 4)


> df <- df[, c('X', 'Y', 'Label')]
> df <- cbind(df,
            data.frame(Gamma1SVM = ifelse(predict(sigmoid.gamma1.svm.fit) > 0, 1, 0),
                       Gamma2SVM = ifelse(predict(sigmoid.gamma2.svm.fit) > 0, 1, 0),
                       Gamma3SVM = ifelse(predict(sigmoid.gamma3.svm.fit) > 0, 1, 0),
                       Gamma4SVM = ifelse(predict(sigmoid.gamma4.svm.fit) > 0, 1, 0)))

> predictions <- melt(df, id.vars = c('X', 'Y'))
> ggplot(predictions, aes(x = X, y = Y, color = factor(value))) +
  geom_point() + facet_grid(variable ~ .)
> ggsave(filename="plot06.png")

f:id:unageanu:20160121150802p:plain

ガンマを変えると決定境界の形が変わります。(としか、書いてくれてない・・・)

アルゴリズムを比較する

SVNと、これまでに学習した、ロジスティック回帰/k近傍法を、同じデータセットに当てはめて比較してみます。

# 対象データの読み込みとクリーニング
> load(file.path('data', 'dtm.RData'))

> set.seed(1)

# 訓練データとテストデータに分割
> training.indices <- sort(sample(1:nrow(dtm), round(0.5 * nrow(dtm))))
> test.indices <- which(! 1:nrow(dtm) %in% training.indices)

> train.x <- dtm[training.indices, 3:ncol(dtm)]
> train.y <- dtm[training.indices, 1]

> test.x <- dtm[test.indices, 3:ncol(dtm)]
> test.y <- dtm[test.indices, 1]

> rm(dtm)

まずは、ロジスティック回帰。

> library('glmnet')
> regularized.logit.fit <- glmnet(train.x, train.y, family = c('binomial'))

最適な Lambda の値を、6章 と同様の手順で求めます。

> lambdas <- regularized.logit.fit$lambda
> performance <- data.frame()

> for (lambda in lambdas) {
  predictions <- predict(regularized.logit.fit, test.x, s = lambda)
  predictions <- as.numeric(predictions > 0)
  mse <- mean(predictions != test.y)
  performance <- rbind(performance, data.frame(Lambda = lambda, MSE = mse))
}

> ggplot(performance, aes(x = Lambda, y = MSE)) +
  geom_point() + scale_x_log10()
> ggsave(filename="lambda.png")

f:id:unageanu:20160121150756p:plain

1e-03 ぐらいがよさそう。 min を使ってMSEが最小になる Lambda を計算します。

# MSEが最小になる `Lambda` を計算。
# 今回のデータだと最小になるものが2つ存在するので、max を使って大きい(=より正則化が厳しい)方を選択している
> best.lambda <- with(performance, max(Lambda[which(MSE == min(MSE))]))

Lambda が決まったので、最終的なMSEを計算します。

> mse <- with(subset(performance, Lambda == best.lambda), MSE)
[1] 0.06830769

0.06 になりました。

次は、線形カーネルSVM

> library('e1071')
> linear.svm.fit <- svm(train.x, train.y, kernel = 'linear')

> predictions <- predict(linear.svm.fit, test.x)
> predictions <- as.numeric(predictions > 0)
> 
> mse <- mean(predictions != test.y)
> mse
[1] 0.128

0.12 で、ロジスティックス回帰よりも悪い結果になりました。

次、ガウスカーネルSVM

> linear.svm.fit <- svm(train.x, train.y, kernel = 'radial')

> predictions <- predict(linear.svm.fit, test.x)
> predictions <- as.numeric(predictions > 0)
 
> mse <- mean(predictions != test.y)
> mse
[1] 0.1421538

今回のデータセットでは、ロジスティック回帰や、線形カーネルSVMよりも悪い結果になりました。 これは、このデータでは理想的な決定境界が線形に近い可能性を示唆しています。

最後は、k近傍法です。

> library('class')

> knn.fit <- knn(train.x, test.x, train.y, k = 50)

> predictions <- as.numeric(as.character(knn.fit))

> mse <- mean(predictions != test.y)
> mse
[1] 0.1396923

kの値を5~50で変えて、一番良い値を使うようチューニングしてみます。

> performance <- data.frame()

# 5~50でkの値を変えて試行
> for (k in seq(5, 50, by = 5)) {
  knn.fit <- knn(train.x, test.x, train.y, k = k)
  predictions <- as.numeric(as.character(knn.fit))
  mse <- mean(predictions != test.y)
  performance <- rbind(performance, data.frame(K = k, MSE = mse))
}

# もっと良いkの値を取り出す。
> best.k <- with(performance, K[which(MSE == min(MSE))])

> best.mse <- with(subset(performance, K == best.k), MSE)
> best.mse
[1] 0.09169231

チューニングの結果、0.09まで改善しました。

ということで、この問題に対してはロジスティック回帰が一番適しているという結論になりました。

以上から、得られる教訓。

  • 実際のデータセットに取り組むときは、複数アルゴリズムを試した方が良い。
  • 最適なアルゴリズムは、問題の構造に依存する。
  • モデルの性能は、パラメータにも依存する。良い結果を得たければパラメータの調整にも時間をかける

感想

  • 回帰分析とか、最適化について、入門くらいはできたかな。
  • 理論もいいけど、具体的にデータを触って、動かして学びたいプログラマ向けには悪くない本かと。(裏表紙にもそのようなコンセプトの本ですと書かれていますし)
  • ただ、Amazonの書評の通り、理論的説明はあまりないので、そこは別の本が必要な感じです。
    • この本の知識だけで Wikipedia とかに行くと ファッ てなります。
  • 機械学習自体は面白かった。これを使って何か作りたい。

おしまい。