読者です 読者をやめる 読者になる 読者になる

六本木で働くデータサイエンティストのブログ

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

Twitterがリリースした時系列異常値検出のためのRパッケージ{AnomalyDetection}を試してみる

もう松の内も明けてしまいましたが、遅ればせながら皆さん明けましておめでとうございます。今年もよろしくお願いいたします。


で、年明け早々にTwitterエンジニアブログに面白いネタが上がっていたのでした。



その名も{AnomalyDetection}というRパッケージ。バルス砲に備えてTwitterが鉄壁の防御を敷いていることは多くの人がご存知だと思いますが(笑)、そういうバルス砲のような異常なアクセスの跳ね上がりだったり、逆にサーバダウンなどでアクセスの異常な落ち込みが出た時にいち早く検出するというのが目的の分析フレームワークのようです。


ということで、話題が新鮮なうちにちょっと試してみようと思います。


どういう仕組みで動いているのか


理論的背景としては、そもそもgeneralized ESD (extremely Studentized deviate) testという外れ値検出の手法があり、これに季節調整とトレンドを加味したというもののようです。なおgeneralized ESD testについてはこちらに資料があります。


詳細はリンク先を見ていただきたいのですが、やり方としては逐次統計量を出しては比較*1という感じのようです。なおREADME.mdの説明には以下のように書いてあります。

How the package works


The underlying algorithm – referred to as Seasonal Hybrid ESD (S-H-ESD) builds upon the Generalized ESD test for detecting anomalies. Note that S-H-ESD can be used to detect both global as well as local anomalies. This is achieved by employing time series decomposition and using robust statistical metrics, viz., median together with ESD. In addition, for long time series (say, 6 months of minutely data), the algorithm employs piecewise approximation - this is rooted to the fact that trend extraction in the presence of anomalies in non-trivial - for anomaly detection.


Besides time series, the package can also be used to detect anomalies in a vector of numerical values. We have found this very useful as many times the corresponding timestamps are not available. The package provides rich visualization support. The user can specify the direction of anomalies, the window of interest (such as last day, last hour), enable/disable piecewise approximation; additionally, the x- and y-axis are annotated in a way to assist visual data analysis.


ちなみにエンジニアブログの方にちょっと書いてありますが、このパッケージの秀逸なところは「大局的に見た場合の異常値も局所的な異常値も逃さない」ところかと。つまり、週次の季節成分から逸脱したようなちょろっとしたスパイクも捉えられるし、もちろんバルス砲みたいなcrazyな巨大スパイクも捉える。これを実現するために季節調整とトレンドを入れたと見るのが妥当そうですね。


なお、本来ならこういう異常値検出フレームワークバッチ処理よりはオンライン処理で動いて欲しい気がするんですが、generalized ESD testは逐次比較型の手法なのでひょっとしたらオンラインでも動くかな?と思いました。まぁ、毎時間ごとにバッチで回しても良いのかも。。。


インストールするには


最近リリースされている新興パッケージ群と同じように、{devtools}を使ってGitHubからコンパイルしつつ入れる感じです。

> install.packages("devtools")
> devtools::install_github("twitter/AnomalyDetection")
> library(AnomalyDetection)


{devtools}が入っている場合は2行目から実行すればOKです。なお、このパッケージの関数で扱えるデータは基本的にはPOSIXltクラスによるタイムスタンプが付与されたデータフレームか、事前に季節周期が判明しているベクトルかのいずれかの模様です。


実際に{AnomalyDetection}をヘルプの例に従って試してみる


README.mdに載ってる例で恐縮ですが、実際にやってみます。これだけで基本的には例と同じ結果が得られるはずです。

> data(raw_data)
> plot(raw_data)
> plot(raw_data,type='l')
> res <- AnomalyDetectionTs(raw_data, max_anoms=0.02, direction='both', plot=TRUE)
> res
$anoms
              timestamp    anoms
1   1980-09-25 16:05:00  21.3510
2   1980-09-29 06:40:00 193.1036
3   1980-09-29 21:44:00 148.1740
4   1980-09-30 17:46:00  52.7478
5   1980-09-30 17:48:00  49.6582
# 中略
127 1980-10-05 01:32:00  47.9571
128 1980-10-05 13:08:00 210.0000
129 1980-10-05 13:18:00  40.0000
130 1980-10-05 13:28:00 250.0000
131 1980-10-05 13:38:00  40.0000

$plot


AnomalyDetectionTs関数の返値は$anomsが検出された異常値、$plotはズバリそのままプロットしてくれるメソッドです。なのでresとやれば異常値を列挙した上にプロットしてくれますし、ヘルプの例だとres$plotとやればプロットだけしてくれます(その代わりplot関数では読み込めない仕様になっている)。プロットの結果は以下の通り。


f:id:TJO:20150109114021p:plain


ということで、うまく出来たっぽいです。


自前のシミュレーションデータでやってみる


ヘルプの例でうまくいくのは当たり前なので、ちょっと色々小細工したデータでもう一度やってみます。これは以前何度か時系列分析まわりで使ったことのある「周期的スパイク+ベースライン変化を伴う」時系列シミュレーションデータです。

> x1<-rnorm(300) # ホワイトノイズ
> x2<-rep(c(15,rep(0,49)),6) # スパイク異常値が6回生じる時系列で周期は50期
> x3<-c(rep(0,150),rep(3,150)) # ステップ状の不連続な変化
> x<-x1+x2+x3
> plot(x)
> y<-x+c(rep(0,130),30,rep(0,169)) # 131期目にサイズ30の異常値を混入させてみた
> res2<-AnomalyDetectionVec(y,max_anoms=0.02,direction='both',plot=TRUE,period=50)
> res2
$anoms
  index    anoms
1   131 29.92602
# ちゃんと131期目の異常値を検出した

$plot

f:id:TJO:20150109114218p:plain


これもうまくいきました。ではちょっと意地悪してみましょう。

> noisex<-rnorm(n=300,mean=0.5,sd=0.2)
> noisex[which(noisex>0.9)]<-5
> noisex[which(noisex<=0.9)]<-0
> y2<-y+noisex
> res3<-AnomalyDetectionVec(y2,max_anoms = 0.02,direction = 'both',period = 50,plot = TRUE)
> res3
$anoms
  index    anoms
1   131 29.92602
# お、サイズ5のスパイクは捕まえられないっぽい

$plot


ダメでしたorz なので、スパイクのサイズを2倍にしてみます。

> noisex<-rnorm(n=300,mean=0.5,sd=0.2)
> noisex[which(noisex>0.9)]<-10
> noisex[which(noisex<=0.9)]<-0
> y3<-y+noisex
> res4<-AnomalyDetectionVec(y3,max_anoms = 0.02,direction = 'both',period = 50,plot = TRUE)
> res4
$anoms
  index      anoms
1    31 -0.2261887
2    45  9.9976652
3   131 29.9260198
4   181  1.2461669
5   281  2.0063777
6   296 12.3916707
# 本来スパイクを入れたはずの点がいくつも見逃され、
# さらに関係ない点が検出されているケースも目立つ

$plot

f:id:TJO:20150109115352p:plain


意外とうまくいってない感がorz という点から察するに、もうちょっと季節性が明確な時系列データでないといけないということなんですかね。。。チャンスがあればちょっと実データにも使ってみようかなと思います。

*1:Johansen's procedureを思い出しますね