はじめに

graphical lasso を用いた異常検知を行ってみた。参考資料は末尾に記載した。今回は時系列、多変量、かつ異常と正常のラベルデータが付随しているデータセットであるSECOM Data Setを使用した。この記事のシリーズでは、

R で graphical lasso を用いた異常検知　その１：基礎分析
R で graphical lasso を用いた異常検知　その２：前処理
R で graphical lasso を用いた異常検知　その３：モデル構築
R で graphical lasso で用いた異常検知　その４：異常度の算出　←　今ここ！

の順番で書いていく。本記事ではモデル構築の結果をもとに、異常解析、外れ値解析を行っていく。以下、コードは 1., 2., 3. の記事のコードが実行済みであるとする。

異常解析の場合の異常度の定義

正常なデータ集合$D$とテストのデータ集合$D'$の間で、どの変数にどれくらい食い違いがあるのかを示すのが異常解析の場合の異常度である。その場合、異常度は以下の様に表される。graphical lasso で求まる精度行列を$\Lambda$とすると、

a_i = \frac{1}{2} \ln \frac{\Lambda_{i, i}}{\Lambda'_{i, i}} - \frac{1}{2} \left\{ \frac{\left[ \Lambda S \Lambda \right]_{i,i}}{\Lambda_{i,i}} - \frac{\left[ \Lambda' S \Lambda' \right]_{i,i}}{\Lambda'_{i,i}} \right\}

となる。$[\cdot ]_{i,i}$は括弧の中の行列の$i,i$成分を取り出す演算を表している。また、$S$は$D$の標本共分散行列である。

以上を用いて、正常データセットと異常データセットの比較を行う。

異常解析結果

glasso で推定したモデルから異常度を推定する関数 abn を定義する。

コード

abn <- function(train, test){
  wi <- train$wi
  w <- train$w
  wi_ <- test$wi
  for(i in 1:ncol(wi)){
    a <- 0.5 * log(wi[i, i]/wi_[i, i]) - 0.5 * ((wi*w*wi)[i, i]/wi[i, i] - (wi_*w*wi_)[i, i]/wi_[i, i])
    return(a)
  }
}

試しに正常データで正常データを推定すると、

コード

normal_test <- abn(fit.norm, fit.norm)
barplot(normal_test, ylab="Abnormality", xlab="Columns")

全て0の結果が表示される。

正常データで異常データのテストを行うと、

コード

abn_test <- abn(fit.norm, fit.abn)
barplot(abn_test, ylab="Abnormality", xlab="Columns")

いくつかの項目で異常度が大きくなる。異常度が大きい項目が異常を示している事が示唆される。いくつか異常度がマイナスを示す項目があるが、そこの解釈は分からなかった。（もし何か知見がある方がいらしたら是非コメントお待ちしております。）

外れ値解析の異常度の定義

正常のデータ（$D$）を使って確率分布 $p(x|D)$ が得られたとする。新しくデータ$x'$が手に入ったときに、$x'$の各変数の異常度は

a_i \left( x \right) = - \ln p \left( x_i' | x'_{-i}, D \right) = \frac{1}{2} \ln \frac{2\pi}{\Lambda_{i,i}} + \frac{1}{2\Lambda_{i,i}}\left( \sum^M_{j=1} \Lambda_{i,j} x'_j \right)^2

となる。

#外れ値解析の結果
疑似的に正常データからモデル構築をし、それを含めた各データ点について外れ値解析を行ってみる。

out <- function(df, l) {
  x <- df %>% select(-label) %>% select(-datetime)
  a <- x
  for (n in 1:nrow(x)) {
    for (i in 1:ncol(x)) {
      u <- 0
      for (j in 1:ncol(x)) {
        u <- u + l[i,j]*x[n,j]
      }
      a[n,i] <- 1/2*log(2*pi/l[i,i]) + (1/(2*l[i,i]))*u^2
    }
  }
  return(a)
}

outliers <- out(d, fit.norm$wi)
outliers_lt <- cbind(data.frame(d$label), data.frame(d$datetime), outliers)
colnames(outliers_lt)[1:2] = c("label", "datetime")
outliers_lt$label <- as.factor(outliers_lt$label)

結果を可視化していく。

dir <- "./plot_out/"
for (i in 3:ncol(outliers_lt)) {
  xn <- as.name(colnames(outliers_lt)[2])
  yn <- as.name(colnames(outliers_lt)[i])
  cl <- as.name("label")
  fname <- paste(dir, colnames(outliers_lt)[i], "_out.png", sep = "")
  
  g <- ggplot(data = outliers_lt, mapping = aes_(x = xn, y = yn, color = cl))
  g <- g + geom_point()
  
  ggsave(file = fname, plot = g, width = 9.44, height = 4, dpi = 300)
}

プロットの例を示すと、

となる。-1が正常データで1が異常データだが、このプロットを見ると異常度によって明確に-1と1が分かれているわけではない。またベースラインが0を示しておらず、オフセットがのっている。この点の解釈についても、もう少し勉強が必要である。

一応、正常、異常の分類問題を解いたことになるので、ROC曲線も描いてみる。

library(pROC)

dir <- "./plot_roc/"
#d <- df.comp.scale.abn.state

auc <- c()
for (i in 3:ncol(outliers_lt)) {
  fname <- paste(dir, colnames(outliers_lt)[i], "_ROC.png", sep = "")
  eval(parse(text=paste0("ROC <- roc(label ~ ", colnames(outliers_lt)[i], ", data = outliers_lt, ci = TRUE)")))
  
  auc <- append(auc, ROC$auc[1])
  
  g <- ggroc(ROC)
  g <- g + ggtitle(paste(colnames(outliers_lt)[i], ":", ROC$auc[1]))
  ggsave(file = fname, plot = g, width = 4, height = 4, dpi = 300)
}

ROCプロットとAUCの値の1例を示すと、

となる。AUC $\simeq$0.58 だったので、ほとんどランダムに判別した（AUC = 0.5）と変わらない結果となった。他の項目についても同じような値を示していた。

おわりに

AUCの値からも言えるように、今回のデータセットでは、graphical model を用いた異常検知で優位な異常の検出には至らなかった。今後、この手法を正しく用いるためにも、自身の実装に問題があったのか、データセットに手法が合わなかったのか、を検討していく必要がある。もし有意義な検討、検証が行えれば、記事を追記したいし、他の異常検知の手法についても調査、実装を行っていきたいと思う。

参考文献（graphical lasso を用いた異常検知について）