Rで異常検知(2): 正規分布に従うデータからの異常検知（ホテリング理論・MT法）

さて、気紛れから始まったこのシリーズですが。今回は第2章を取り上げます。

作者: 井手剛
出版社/メーカー: コロナ社
発売日: 2015/02/19
メディア: 単行本
この商品を含むブログ (4件) を見る

多変量かつ非正規データの異常検知は少し後の方になるので、例のwater treatment plantのデータセットを持ち出すのは後回しにして、今回は適当に生成したデータセットを使うことにしました。

ちなみに、今回のシリーズではあまりあれこれ引用しまくると引用の範囲を超えてしまいそうな気がしたので（笑）、要点をちろっとまとめてRスクリプトを並べるだけに留めておきます。故に、皆さんご自身がお手元で試される場合には必ず井手先生のテキストをご用意下さい、ということで。

ホテリングの $T^2$ （1次元）

これはオーソドックスに正規分布する1次元のデータの中から、全体の分布に対して逸脱している異常値を検出するタイプの方法論です。

手順自体はテキストのpp.23-25に書いてありますが、井手さんが何度も強調されているように「そこに至るまでの理論的な積み上げ」を知ることの方が重要なので、pp.19-22及びpp.26-36の理論面の解説にも目を通しておくことをお薦めしておきます。

手順2.1
0) 準備
異常が含まれていないか、含まれていたとしてもごく少数と思われるデータセットを用意する。異常判定の閾値を確率値 $\alpha$ で与え（例えば0.01や0.03）、カイ二乗分布の表から、異常度の閾値 $\alpha_{th}$ を求めておく。
1) ステップ1（分布推定）
標本平均および標本分散を計算する。
2) ステップ2（異常度の計算）
新たな観測値 $x'$ が得られるたび、異常度
$a(x') = \left( \frac{x' - u}{\hat{\sigma}}\right)^2$
の値を計算する。
3) ステップ3（閾値判定）
異常度が閾値 $\alpha_{th}$ を超えたら異常と判定する。

ということでこの通りにやってみます。なお、以下で生成したデータセットは以前{AnomalyDetection}パッケージを紹介した記事の中で生成したのと同じものです。

> x1<-rnorm(300) # ホワイトノイズ
> x2<-rep(c(15,rep(0,49)),6) # サイズ15のスパイク異常値が6回生じる時系列で周期は50期
> x3<-c(rep(0,150),rep(3,150)) # ステップ状の不連続な変化
> x<-x1+x2+x3
> y<-x+c(rep(0,130),30,rep(0,169)) # 131期目にサイズ30の異常値を混入させてみた
> plot(y, type='l')

> mu <- mean(y) # 標本平均
> s2 <- mean((y - mu)^2) # 標本分散（不偏分散ではない）
> mu
[1] 1.806073
> s2
[1] 10.27356
> a <- (y - mu)^2 / s2 # 異常度
> thr <- qchisq(0.99, 1) # カイ二乗分布から99%閾値を求める
> plot(a) # 異常度をプロット
> segments(0, thr, 300, thr, col='red', lty=3, lwd=3) # 閾値に線を引く

f:id:TJO:20170208154651p:plain
f:id:TJO:20170205164320p:plain

確かにスパイクが全て閾値を超えていて異常値と判定されています。面白いのが、1個だけ紛れ込ませたサイズ30の異常値の異常度が、その他の周期的に発生するサイズ15のスパイクよりも絶対値で比較した場合よりもより大きく出ているところ。異常度はある程度「より異常かどうか」も反映するということなんでしょうか。

ホテリングの $T^2$ （多次元）

先ほどの方法が正規分布に従う1次元データにおける異常検知であるならば、こちらは多次元データにおける同じ方法論です。

こちらも手順自体はp.42に書いてありますが、pp.37-41の理論の導入及びpp.44-48の理論の詳細についても目を通しておくことをお薦めします。大事なことはマハラノビス距離を異常度として採用している点と、Mが大きい場合の漸近性が仮定されている点なのかな？と個人的には理解しています。

手順2.2
あらかじめある所与のパーセント値 $\alpha$ に基づき、自由度M・スケール因子1のカイ二乗分布から方程式
$1 - \alpha = \displaystyle \int_0^{a_{th}} du ~ \chi^2 (u | M, 1 )$
により閾値 $a_{th}$ を求めておく。
1) 正常標本が圧倒的多数を占めると信じられるデータから標本平均および標本共分散行列
$\hat{\mathbf{\mu}} = \frac{1}{N} \displaystyle \sum^N_{n=1} \mathbf{x}^{(n)}$ ,　 $\hat{\mathbf{\Sigma}} = \frac{1}{N} \displaystyle \sum^N_{n=1} (\mathbf{x}^{(n)} - \hat{\mathbf{\mu}}) (\mathbf{x}^{(n)} - \hat{\mathbf{\mu}})^T$
を計算しておく。
2) 新たな観測値 $\mathbf{x}'$ を得るたび、異常度としてのマハラノビス距離
$a(\mathbf{x}') = (\mathbf{x}' - \hat{\mathbf{\mu}})^T \hat{\mathbf{\Sigma}}^{-1} (\mathbf{x}' - \hat{\mathbf{\mu}})$
を計算する。
3) $a(\mathbf{x}') > a_{th}$ なら警報を発する。

こちらの実践なんですが、テキストと同じように2次元でやっても面白くないので、ギリギリ可視化可能な3次元のデータを3次元正規分布から適当に生成して試してみます。

# 3次元正規分布からデータを生成したつもり
> x <- rnorm(n = 500, mean = 1, sd = 1)
> y <- rnorm(n = 500, mean = 1, sd = 1)
> z <- rnorm(n = 500, mean = 1, sd = 1)
> d <- data.frame(x = x, y = y, z = z)

# 異常値を適当に3つ作って一応可視化しておく
> idx <- sample(nrow(d), 3)
> d[idx,] <- d[idx,] + 10
> install.packages("scatterplot3d")
> library(scatterplot3d)
> scatterplot3d(d, pch=19, cex.symbols=0.5)

# ホテリングのT^2を算出する
> X <- as.matrix(d)
> mx <- colMeans(X)
> Xc <- as.matrix(X) - matrix(1, nrow(X), 1) %*% mx
> Sx <- t(Xc) %*% Xc / nrow(X)
> am <- rowSums((Xc %*% solve(Sx)) * Xc)
> thr <- qchisq(0.99, 3)
> plot(am)
> segments(0, thr, 500, thr, col='red', lty=3, lwd=3)

# 異常値と判定されたデータ点を赤く大きく表示する
> idx_a <- which(am > thr)
> col_3d <- rep(1, 500)
> col_3d[idx_a] <- 2
> cex_3d <- rep(0.5,500)
> cex_3d[idx_a] <- 2
> scatterplot3d(d, pch=19, cex.symbols=cex_3d, color=col_3d)

f:id:TJO:20170424120321p:plain
f:id:TJO:20170424120341p:plain
f:id:TJO:20170424120400p:plain

思ったよりも異常と判定されたデータ点が多くなりましたが、大体こんな感じかなという結果になりました。なお、試しにwater treatment plantのデータセットでやるとこうなります。

> d <- read.csv('watertreatment_mod.csv')
> d1 <- d[,-1]
> mx <- colMeans(d1)
> Xc <- as.matrix(d1) - matrix(1, nrow(d1), 1) %*% mx
> Sx <- t(Xc) %*% Xc / nrow(d1)
> am <- rowSums((Xc %*% solve(Sx)) * Xc)
> thr <- qchisq(0.99, ncol(d1))
> plot(am)
> segments(0, thr, 380, thr, col='red', lty=3, lwd=3)

f:id:TJO:20170424120638p:plain

明らかに異常度のスケールがおかしくなっています。これではどうにもならないわけで、今後紹介される別の方法を使うべきだということなのでしょう。

マハラノビス=タグチ法（MT法）

これまでの方法は単純にあるデータ点（データマトリクスの中では特定の「行」）が異常か否かを判別してきました。これに対し、マハラノビス=タグチ法（MT法）は「異常値の判別に貢献する特徴量はどれか」をS/N比の概念を用いて示すものです。テキストpp.49-52に手順が載っていますが、これまた

手順2.3
正常データが圧倒的多数だと信じられるデータセット $\cal{D} = \{ \mathbf{x}^{(1)}, \cdots, \mathbf{x}^{(N)} \}$ と、異常と判明しているデータセット $\cal{D}' = \{ \mathbf{x}'^{(1)}, \cdots, \mathbf{x}'^{(N)} \}$ を用意する。変数の数は共に $M$ （したがってデータは $M$ 次元ベクトルの集まり）とする。
1) 分布推定
$\cal{D}$ を基に、標本平均と標本共分散を求める。
2) 異常度の計算
$\cal{D}$ の中の各標本に対し、1変数当たりのマハラノビス距離（異常度を $M$ で割ったもの）を計算する。
3) 異常判定1
$\cal{D}$ の標本が正常範囲に入るように1変数当たりのマハラノビス距離の閾値を決める。
4) 異常判定2
$\cal{D}'$ の各標本に対して、 $M$ 変数の中からいくつかの変数を選び、その変数集合の1変数当たりの異常度を計算する。

やっていることは多次元のホテリングの $T^2$ の途中までと共通で、異常度だけ特徴次元数で割っておくというものです。閾値は1としておけば良い、ということになっています。

> # 3次元正規分布からデータを生成したつもり
> xt <- rnorm(n = 500, mean = 1, sd = 1)
> yt <- rnorm(n = 500, mean = 1, sd = 1)
> zt <- rnorm(n = 500, mean = 1, sd = 1)
> dt <- data.frame(x = xt, y = yt, z = zt)

> # 異常値をz軸方向に逸脱するように適当に3つ作って一応可視化しておく
> idxt <- sample(nrow(dt), 3)
> dt[idxt,3] <- dt[idxt,3] + 20
> library(scatterplot3d)
> scatterplot3d(dt, pch=19, cex.symbols=0.5)

> # ホテリングのT^2とMT法による異常度を算出する
> Xt <- as.matrix(dt)
> mxt <- colMeans(Xt)
> Xct <- as.matrix(Xt) - matrix(1, nrow(Xt), 1) %*% mxt
> Sxt <- t(Xct) %*% Xct / nrow(Xt)
> amt <- rowSums((Xct %*% solve(Sxt)) * Xct) / ncol(Xt)
> thrt <- 1
> plot(amt)
> segments(0, thrt, 500, thrt, col='red', lty=3, lwd=3)

> # 異常値と判定されたデータ点を赤く大きく表示する
> idx_at <- which(amt > thrt)
> col_3dt <- rep(1, 500)
> col_3dt[idx_at] <- 2
> cex_3dt <- rep(0.5,500)
> cex_3dt[idx_at] <- 2
> scatterplot3d(dt, pch=19, cex.symbols=cex_3dt, color=col_3dt)

f:id:TJO:20170208170621p:plain
f:id:TJO:20170208170641p:plain
f:id:TJO:20170208170658p:plain

z軸方向に異常値を生成したので、z座標が恐らく最も異常値を良く表現しているのではないかと考えられます。そこで、p.52に従って以下のようにS/N比を計算してみましょう。

> par(mfrow = c(2,2))
> for (i in 1:3){
+     xc_prime <- Xct[idxt[i],]
+     SN1 <- 10*log10(xc_prime^2 / diag(Sxt))
+     barplot(SN1)
+ }

f:id:TJO:20170208171229p:plain

負のS/N比は平均からの偏差が標準偏差より小さいことを意味しているので、「逸脱」という意味での異常検知を行う上では特に関係がありません。その点で言うと、z座標だけが一貫して絶対値の大きい正のS/N比を示しており、このケースではやはり異常度はz座標の値に主に帰せられるということが見て取れます。

最後に

とりあえず井手先生のテキストに沿ってザッとなぞってみました。次回は3章の「非正規データからの異常検知」をやってみようと思いますが、この3章が結構範囲が広くて長いので何回かに分けるつもりです。

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

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

Rで異常検知(2): 正規分布に従うデータからの異常検知（ホテリング理論・MT法）

ホテリングの $T^2$ （1次元）

ホテリングの $T^2$ （多次元）

マハラノビス=タグチ法（MT法）

最後に

ホテリングの（1次元）

ホテリングの（多次元）

マハラノビス=タグチ法（MT法）

最後に

ホテリングの $T^2$ （1次元）

ホテリングの $T^2$ （多次元）