ガウス過程回帰・分類をRで試してみた

先日こちらの書籍をご恵贈いただきました。

ガウス過程と機械学習 (機械学習プロフェッショナルシリーズ)

ガウス過程と機械学習ということで、本書の6.2節でも取り上げられていてハイパーパラメータ探索法として用いられることもあるベイズ最適化を僕は真っ先に思い出したのですが、以下の過去エントリを書いた際に「これって何でこうやってるんだろう」と思っていたことの大半が書かれていて、「なるほど」と納得すると同時に自分の勉強不足を恥じた次第です。。。「ガウス過程回帰法の特長をフル活用した応用技術」だったんですね。

ということで少しずつ読み進めているんですが、そう言えばRでガウス過程回帰出来るパッケージなかったっけ？と思っていたらありました。だいぶ老舗のパッケージでもある{kernlab}と、その他いくつかのCRAN / GitHub配布パッケージにガウス過程回帰・分類のメソッドが実装されているようです。そんなわけで、適当なサンプルデータを用意してサクッとガウス過程回帰・分類を試してみようと思います。ちなみにPythonの場合は普通にscikit-learnにgaussian_processがあるので、そこからimportすれば簡単に実践できます。こちらの方はググれば実践例が色々見つかりますので、今回は割愛します。

そしていつもながらですが、雑に勉強しながら適当に実践してみただけという記事ですので、誤りや理解不足のところなどあればガンガン突っ込んでいただけると有難いですm(_ _)m

ガウス過程回帰とは

(Gaussian process - Wikipedia)
言うまでもなくご恵贈いただいた上記の書籍をお読みいただくのが宣伝にもなってベストなのですが（笑）、時間がないという方のために超絶簡潔な説明を。

これは僕の個人的な理解なのですが、一言で表せば「あらかじめ何かしらのモデル（線形回帰モデルなど）を定めてそれとデータとを比べて最適化してパラメータを求めることでモデルを推定するのではなく、モデル全体をガウシアンカーネルなどの『カーネル』の多数（それこそ無限）の集合体とみなして柔軟にフィッティングさせてモデルとする」手法のようです。これは当然計算コストが高くなりますが、積分消去によって周辺化することである程度下げることができます。太字にしましたが、ガウス過程回帰（そして分類）はその柔軟さが最大の長所であり、そのために重宝される局面が存在するという理解です。

あとは、今回は試しませんでしたがガウス過程回帰はその性質上、予測値の平均と分散を算出することも出来ます。現在のところこれが大きな利点となる局面は時系列予測以外ではちょっと思いつきませんが。。。*1

多項式フィッティングをしてみる

ということで、早速ちょっとした例で試してみましょう。まずは定番の多項式フィッティングですが、これは単純に以下のコードを回すだけです。GitHubにも上げてあります。gausspr関数はtype引数を与えなければ回帰(regression)が走ります。紫で強調した点が学習データで、残りがテストデータです。

また、ただgausspr回すだけでは面白くないので{e1071}のsvmでも同じ多項式フィッティングをして、どう変わるかを比較してあります。

library(kernlab)
x <- seq(-10, 10, length.out = 300)
set.seed(51)
y <- -1e-3 * x * (x + 9) * (x - 1) * (x - 9) + rnorm(300, 0, 0.5)
y_true <- -1e-3 * x * (x + 9) * (x - 1) * (x - 9)
plot(x, y, xlim = c(-11, 11), ylim = c(-3, 3), xlab = '', ylab = '')
lines(x, y_true, type = 'l', col = 'blue', lwd = 2)
d <- data.frame(x = x, y = y)
set.seed(71)
d1 <- d[sample(300, 50, replace = F),]
fit <- gausspr(y~x, d1)
fit
y_pred <- predict(fit, d)
par(new = T)
plot(d1$x, d1$y, pch = 19, col = 'purple',
     xlim = c(-11, 11), ylim = c(-3, 3), xlab = '', ylab = '')
lines(x, y_pred, type = 'l', col = 'red', lwd = 3)

library(e1071)
fit_svm <- svm(y~x, d1)
y_pred_svm <- predict(fit_svm, d)
lines(x, y_pred_svm, type = 'l', col = '#008000', lwd = 3)

legend('bottom', legend = c('True', 'GP', 'SVM'),
       col = c('blue', 'red', '#008000'), lwd = c(2, 3, 3))

> fit <- gausspr(y~x, d1)
Using automatic sigma estimation (sigest) for RBF or laplace kernel 
> fit
Gaussian Processes object of class "gausspr" 
Problem type: regression 

Gaussian Radial Basis kernel function. 
 Hyperparameter : sigma =  4.55859552003361 

Number of training instances learned : 50 
Train error : 0.324707014 

両端がちょっと当たっていない気もしますが、それ以外のところは大体フィットしているようにも見えます。GPの方がSVMよりも僅かに端側で良くフィットしているようです。

XORパターン分類をしてみる

ガウス過程回帰は分類確率を出すように変えれば分類器としてモデリングすることも出来ます。ここではGitHubに上げてあるXORデータセットに対して分類を行ってみます。

多項式フィッティングの時と同様、比較対象として{e1071}のsvmによる分類も同時に行っています。

library(kernel)
d <- read.csv('https://raw.githubusercontent.com/ozt-ca/tjo.hatenablog.samples/master/r_samples/public_lib/jp/xor_complex_small.txt', sep = ' ')
summary(d)
d$label <- as.factor(d$label-1)
fit <- gausspr(label~., d, type = 'classification')
out <- predict(fit, newdata = pgrid, type = 'probabilities')
out_vec <- rep(0, nrow(out))
for (i in 1:nrow(out)){
  if (out[i, 1] > 0.5) out_vec[i] <- 1
}
plot(c(), type = 'n', xlim = c(-4, 4), ylim = c(-4, 4), xlab = '', ylab = '')
par(new = T)
rect(0, 0, 4, 4, col = '#dddddd')
par(new = T)
rect(-4, 0, 0, 4, col = '#ffdddd')
par(new = T)
rect(-4, -4, 0, 0, col = '#dddddd')
par(new = T)
rect(0, -4, 4, 0, col = '#ffdddd')
par(new = T)
plot(d[, -3], pch = 19, cex = 2, col = d$label, xlim = c(-4, 4), ylim = c(-4, 4),
     xlab = '', ylab = '')
par(new = T)
contour(px, py, array(out_vec, c(length(px), length(py))),
        col = 'purple', levels = 0.5, lwd = 5, drawlabels = F)

library(e1071)
fit_svm <- svm(label~., d)
out_svm <- predict(fit_svm, newdata = pgrid)
par(new = T)
contour(px, py, array(out_svm, c(length(px), length(py))),
        col = '#008000', levels = 0.5, lwd = 5, drawlabels = F)
legend('topleft', legend = c('GP', 'SVM'), col = c('purple', '#008000'), lwd = 5)

それっぽい感じになりました。サンプルサイズが100だと、GP / SVMのどちらとも相応に歪んだ感じの決定境界を示します。

ちなみにサンプルサイズを1万まで増やしたバージョン（GitHubで "~medium.txt" という名前になっているやつです）に置き換えるとこうなります。GPとSVMとでほぼ同じ決定境界を描いているように見えます。が、GPの方の計算時間はかなり長かったです。。。

感想

実は既にこういうQiita記事がだいぶ前に出ていまして。

感想としては大体一緒です。XORパターン分類で1万サンプルサイズのデータに対してGPだとまともな時間では全然回り切らず、一方でSVMは（以前試した感じでは）それほど莫大な時間もかからずに回り切っていたので、計算コストという点ではSVMの方が有利なように見えます。

これに対して、過学習しにくいというか汎化性能の高さという点では回帰においてはGPの方が良いかも？と思いましたが、分類においては正直微妙な感じがしました。ただ、分類ではSVMが「丸まっていく」感じの決定境界を描いている（なので図の外では円弧を描く）一方でGPは多少ぐにゃぐにゃしながらも直角になる方向に向かっているようにも見え、何とも言えない気もします。

ということで、もうちょっとGPの有難みが良く分かる課題で改めて試してみようかなと思った次第です。。。お粗末様でした。

*1:そして予測値の信頼区間が欲しいだけなら従来の自己回帰系の手法でも十分だったりする