（注：ただの備忘録ゆえ、ほぼ確実に後で追記が出る見込みです）
今回はコロナ社井手本の第3章を取り上げます。
- 作者: 井手剛
- 出版社/メーカー: コロナ社
- 発売日: 2015/02/19
- メディア: 単行本
- この商品を含むブログ (4件) を見る
とは言っても全部丸写しするのはさすがに問題があるので、個人的に興味のあるトピックスだけを取り上げて、なおかつあくまでも個人的な備忘録として記する程度に留めておきます。本格的に勉強したい方はちゃんとコロナ社井手本をお求めの上ご自身で独学するようにして下さい。
ということで、今回個人的に取り上げるのは第3章の中のカーネル密度推定とクラスタリングとSVMの箇所です。なおコロナ社井手本ではここでEMアルゴリズムそのものの説明もなされていますが、EMアルゴリズム自体はこのブログでも一度取り上げているので既知のものとしてこの記事では取り扱います。
ただし、個人的にはコロナ社井手本p.71に載っている混合モデル＋EMアルゴリズムの簡単な実践例は割と親切で良いというか、多分知っている限りで最も分かりやすいEMアルゴリズムの解説だと思うので、EMアルゴリズム自体が初見という人はそこの前後の箇所をよく読んで勉強してみるのをお薦めいたします。少なくとも僕には非常に良い復習になりました。
ちなみに今回の記事はちょっと色々手が回らなかったのもあって結構手抜き気味です。後で追記の形で勉強し直した結果などが加わる可能性がありますので、どうしても気になるという方は後ほど改めてご覧ください。。。
カーネル密度推定（※異常検知はできなかった）
井手本pp.78-83の内容です。ここでは{KernSmooth}パッケージを使って、テキストとは異なりド定番のOld Faithful Geyserデータセットを対象に試してみようと思います。ただし、p.83の異常度の評価のところのカーネル行列Kがどこで出てきたのか（もしくはどう入れるべきか）なのかが分からず。。。一旦ここではスキップしてあります。やり方が分かったら追記します。
> data("faithful") > d <- faithful > library(KernSmooth) > h <- c(dpik(d$eruptions), dpik(d$waiting)) > est <- bkde2D(d, bandwidth = h, gridsize = c(500, 500)) > y <- list(x=est$x1, y=est$x2, z=est$fhat) > image(y, col=terrain.colors(7)) > contour(y, add=T)
大体元のcarデータセットでやるのと似たような感じになりました。
クラスタリング（混合正規分布＋EMアルゴリズム）
これは完全に井手本をそのままなぞる形で、データセットをOld Faithful Geyserに変えただけです。
> library(mclust) > res <- Mclust(d) > print(summary(res, parameters = TRUE)) ---------------------------------------------------- Gaussian finite mixture model fitted by EM algorithm ---------------------------------------------------- Mclust EEE (ellipsoidal, equal volume, shape and orientation) model with 3 components: log.likelihood n df BIC ICL -1126.361 272 11 -2314.386 -2360.865 Clustering table: 1 2 3 130 97 45 Mixing probabilities: 1 2 3 0.4632682 0.3564512 0.1802806 Means: [,1] [,2] [,3] eruptions 4.475059 2.037798 3.817687 waiting 80.890383 54.493272 77.650757 Variances: [,,1] eruptions waiting eruptions 0.07734049 0.4757779 waiting 0.47577787 33.7403885 [,,2] eruptions waiting eruptions 0.07734049 0.4757779 waiting 0.47577787 33.7403885 [,,3] eruptions waiting eruptions 0.07734049 0.4757779 waiting 0.47577787 33.7403885 > plot(res) Model-based clustering plots: 1: BIC 2: classification 3: uncertainty 4: density Selection: 1 Model-based clustering plots: 1: BIC 2: classification 3: uncertainty 4: density Selection: 2 Model-based clustering plots: 1: BIC 2: classification 3: uncertainty 4: density Selection: 3 Model-based clustering plots: 1: BIC 2: classification 3: uncertainty 4: density Selection: 4 Model-based clustering plots: 1: BIC 2: classification 3: uncertainty 4: density Selection: 0
> pi <- res$parameters$pro > dd <- cdens(modelName = res$modelName, d, parameters = res$parameters) > a <- -log(as.matrix(dd) %*% as.matrix(pi)) > plot(a, ylab='anomaly score')
うぎゃー、これだと全然異常検知できてないですね。ってか、これはそもそも「外れ値なんてこのデータセットにはない」という話になるんでしょうか。。。ということで諦めてwater treatment plantデータセットでやってみます。
> d1 <- read.csv('watertreatment_mod.csv') > d1 <- d1[,-1] > res2 <- Mclust(d1) # 普通のマシンだといつまで経っても終わらないかもしれない > pi2 <- res2$parameters$pro > d2 <- cdens(modelName = res2$modelName, d1, parameters = res2$parameters) > a2 <- -log(as.matrix(d2) %*% as.matrix(pi2)) > plot(a2, ylab='anomaly score' > head(sort(a2, decreasing = T), n=3) [1] 273.3771 208.1646 197.4208 > which(a2 > 197) [1] 4 64 301
試しに上位3つを取ってみて、4, 64, 301のインデックスのサンプルが異常値っぽいということが分かりました。
K-meansで我流で異常検知
これはこのシリーズ記事の第1回でお見せしたのと全く同じやり方をしただけです。K-meansで徐々に想定クラスタ数を増やしていって、ある一定ラインから同じ個数のサンプルしか含まないクラスタが続くようになったらそこで打ち切って、その中で孤立したクラスタを取ってくるという我流のやり方です。同じくwater treatment plantデータセットに対してやってみます。
> d1 <- read.csv('watertreatment_mod.csv') > for (i in 4:10){ + km <- kmeans(d1[,-1], centers=i) + print(table(km$cluster)) + } 1 2 3 4 104 176 39 61 1 2 3 4 5 75 144 50 17 94 1 2 3 4 5 6 16 116 82 67 47 52 1 2 3 4 5 6 7 67 3 52 81 49 112 16 1 2 3 4 5 6 7 8 84 3 64 16 40 43 67 63 1 2 3 4 5 6 7 8 9 30 24 32 47 48 13 81 102 3 1 2 3 4 5 6 7 8 9 10 43 74 44 30 56 24 61 13 32 3
> for (i in 8:10){ + km <- kmeans(d1[,-1], centers=i) + cls <- which(table(km$cluster)==3) + print(which(km$cluster==cls)) + } [1] 139 149 371 [1] 139 149 371 [1] 139 149 371
サンプル3個が分類され続けるクラスタがで見つかったので、それをチェックしてみたところ。。。あれれ、EMアルゴリズムでやった結果と一つも一致しませんね。これはどういうことなんだろう。。。多分僕が間違っているんじゃないかと思いますので、解決したらそのうち追記します（汗）。
1クラスSVM
個人的には今回ちゃんとやりたかったのがこの1クラスSVM。時系列の異常検知はこれまでにも何度かこのブログでも話題にしてきましたが、普通のクロスセクションデータに対する異常検知はあまり真面目にやったことがなく、一方で実務で1クラスSVMを使っているという話を外で割と耳にするので、気になっていたのでした。
で、井手本では{kernlab}でやっていますが、これをそのままなぞってもただコピーしてきただけになりますので（汗）、個人的に愛用している{e1071}で同様のことをやってみようと思います。ただし、ここで使っているwater treatment plantデータセットは既に異常値が混ざったものなのでこのままでは使えません。そこで、上記のEMアルゴリズムで検出された異常値インデックスの4, 64, 301とK-meansで検出されたインデックス139, 149, 371の6サンプルを異常値とみなし、残りを正常値とします。そして正常値データセット「のみ」で1クラスSVMを学習させて正常値の範囲となる超平面を決定させ、それを用いて異常値を判定させてみてどうなるかを見ます。
なお、間違っていると嫌だなぁと思ったので、はじパタpp.131-133でも確認しておきました。個人的には2クラス（多クラス）SVMは散々やってきた上にスクラッチから組んだこともあるので大体知っていたのですが、1クラスSVMの理屈はよく分かっていなかったので勉強になりました。
- 作者: 平井有三
- 出版社/メーカー: 森北出版
- 発売日: 2012/07/31
- メディア: 単行本（ソフトカバー）
- 購入: 1人 クリック: 7回
- この商品を含むブログ (3件) を見る
ちなみに今回試してみて知ったのですが、svm {e1071}の1クラスSVMは異常値だとFALSEを返すんですね。。。一瞬ドキッとしました*1。引数type = 'one-classification'とすれば残りのパラメータ設定は通常の2クラス（多クラス）のC-SVCと同じで大丈夫です。
> train <- d1[-c(4, 64, 139, 149, 301, 371), -1] > test <- d1[c(4, 64, 139, 149, 301, 371), -1] > label <- rep(TRUE, nrow(train)) > tuned <- tune.svm(train, y=label, type='one-classification', nu=0.01:0.99) > tuned$best.parameters nu 1 0.01 > train.one.svm <- svm(train, y=NULL, type='one-classification', nu=0.01) > pred <- predict(train.one.svm, test) > table(pred) pred FALSE 6
とりあえず6個とも異常値と判定されました。しかし、上記の通りEMアルゴリズムとK-meansとでは異なる異常値インデックスを返しています。そこで、それぞれ別々に1クラスSVMで判定してみます。
> # EMアルゴリズムの結果に基づく > train <- d1[-c(4, 64, 301), -1] > test <- d1[c(4, 64, 301), -1] > label <- rep(TRUE, nrow(train)) > tuned <- tune.svm(train, y=label, type='one-classification', nu=0.01:0.99) > tuned$best.parameters nu 1 0.01 > train.one.svm <- svm(train, y=NULL, type='one-classification', nu=tuned$best.parameters$nu) > pred <- predict(train.one.svm, test) > table(pred) pred FALSE 3 > # K-meansの結果に基づく > train <- d1[-c(139, 149, 371), -1] > test <- d1[c(139, 149, 371), -1] > label <- rep(TRUE, nrow(train)) > tuned <- tune.svm(train, y=label, type='one-classification', nu=0.01:0.99) > tuned$best.parameters nu 1 0.01 > train.one.svm <- svm(train, y=NULL, type='one-classification', nu=tuned$best.parameters$nu) > pred <- predict(train.one.svm, test) > table(pred) pred FALSE 3
ということで、とりあえずどちらに基づいた1クラスSVMモデルであってもきちんと異常値判定できました。ただ、これには問題がありまして。。。試しにtrain（学習データ）自身について異常値判定させてみると、こうなります。
> train <- d1[-c(4, 64, 139, 149, 301, 371), -1] > test <- d1[c(4, 64, 139, 149, 301, 371), -1] > label <- rep(TRUE, nrow(train)) > tuned <- tune.svm(train, y=label, type='one-classification', nu=0.01:0.99) > tuned$best.parameters nu 1 0.01 > train.one.svm <- svm(train, y=NULL, type='one-classification', nu=tuned$best.parameters$nu) > pred <- predict(train.one.svm, train) > table(pred) pred FALSE TRUE 32 342
異常値が32個も出てきてしまいました。ちなみにカーネルトリックを使っているということを考えれば当たり前なんですが*2、nuを小さくしたからと言って異常値が減るわけではなくて、例えばnu = 0.0001にすると異常値は130個にまで増えてしまいます。そこで、ガウシアンカーネルのパラメータであるgammaもチューニングで調整してみることにします。
> train <- d1[-c(4, 64, 139, 149, 301, 371), -1] > test <- d1[c(4, 64, 139, 149, 301, 371), -1] > label <- rep(TRUE, nrow(train)) > tuned <- tune.svm(train, y=label, type='one-classification', nu=0.01:0.99, gamma=1e-4:1) > tuned$best.parameters gamma nu 1 1e-04 0.01 > train.one.svm <- svm(train, y=NULL, type='one-classification', nu=tuned$best.parameters$nu, gamma=tuned$best.parameters$gamma) > pred <- predict(train.one.svm, train) > table(pred) pred FALSE TRUE 5 369
5個が異常値と判断されたようで、何となくそれっぽい結果になりました。が、よーく調べてみると、
> table(predict(train.one.svm, test)) FALSE TRUE 2 4
あれ、testの6個中2個しか異常値と判定されていないorz なお、これもgamma = 1e-7とかすると異常値が爆増して訳が分からないことになります。この辺の挙動を理解するにはもっと実装・手法の調査と理論の勉強が必要そうです。。。
反省
全体として見ると今回はうまくいかないところだらけで、完全にただの備忘録になりました。。。もっと色々勉強します。
