渋谷駅前で働くデータサイエンティストのブログ

元祖「六本木で働くデータサイエンティスト」です / 道玄坂→銀座→東京→六本木→渋谷駅前

Rで異常検知(3): 非正規データからの異常検知(カーネル密度推定・EMアルゴリズム・K-means・1クラスSVM)

(注:ただの備忘録ゆえ、ほぼ確実に後で追記が出る見込みです)

今回はコロナ社井手本の第3章を取り上げます。

とは言っても全部丸写しするのはさすがに問題があるので、個人的に興味のあるトピックスだけを取り上げて、なおかつあくまでも個人的な備忘録として記する程度に留めておきます。本格的に勉強したい方はちゃんとコロナ社井手本をお求めの上ご自身で独学するようにして下さい。


ということで、今回個人的に取り上げるのは第3章の中のカーネル密度推定とクラスタリングSVMの箇所です。なおコロナ社井手本ではここでEMアルゴリズムそのものの説明もなされていますが、EMアルゴリズム自体はこのブログでも一度取り上げているので既知のものとしてこの記事では取り扱います。

ただし、個人的にはコロナ社井手本p.71に載っている混合モデル+EMアルゴリズムの簡単な実践例は割と親切で良いというか、多分知っている限りで最も分かりやすいEMアルゴリズムの解説だと思うので、EMアルゴリズム自体が初見という人はそこの前後の箇所をよく読んで勉強してみるのをお薦めいたします。少なくとも僕には非常に良い復習になりました。


ちなみに今回の記事はちょっと色々手が回らなかったのもあって結構手抜き気味です。後で追記の形で勉強し直した結果などが加わる可能性がありますので、どうしても気になるという方は後ほど改めてご覧ください。。。


カーネル密度推定(※異常検知はできなかった)


井手本pp.78-83の内容です。ここでは{KernSmooth}パッケージを使って、テキストとは異なりド定番のOld Faithful Geyserデータセットを対象に試してみようと思います。ただし、p.83の異常度の評価のところのカーネル行列Kがどこで出てきたのか(もしくはどう入れるべきか)なのかが分からず。。。一旦ここではスキップしてあります。やり方が分かったら追記します。

手順3.4(カーネル密度推定におけるバンド幅選択)

  1.  Hの候補値をあらかじめH^1, H^2, \cdotsと用意する
  2.  Hの候補それぞれについて積分二乗誤差(3.35)の値を計算し、記録する
  3. 最小の積分二乗誤差を与える H^iを最適解として採用する
> 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)

f:id:TJO:20170519123303p:plain

大体元の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

f:id:TJO:20170519123606p:plain

f:id:TJO:20170519123619p:plain

f:id:TJO:20170519123647p:plain

f:id:TJO:20170519123712p:plain

> 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')

f:id:TJO:20170519124129p:plain

うぎゃー、これだと全然異常検知できてないですね。ってか、これはそもそも「外れ値なんてこのデータセットにはない」という話になるんでしょうか。。。ということで諦めて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

f:id:TJO:20170519135239p:plain


試しに上位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個が分類され続けるクラスタ k \geq 8で見つかったので、それをチェックしてみたところ。。。あれれ、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の理屈はよく分かっていなかったので勉強になりました。

はじめてのパターン認識

はじめてのパターン認識

ちなみに今回試してみて知ったのですが、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とかすると異常値が爆増して訳が分からないことになります。この辺の挙動を理解するにはもっと実装・手法の調査と理論の勉強が必要そうです。。。


反省


全体として見ると今回はうまくいかないところだらけで、完全にただの備忘録になりました。。。もっと色々勉強します。

*1:試しに全説明変数の平均をcolMeansしたサンプルを作って分類するとTRUEが返るのでそれで分かる

*2:もちろんradialのままなのかpolynomialの方がいいのかとか色々迷いどころは他にもある