もう松の内も明けてしまいましたが、遅ればせながら皆さん明けましておめでとうございます。今年もよろしくお願いいたします。
で、年明け早々に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関数では読み込めない仕様になっている)。プロットの結果は以下の通り。
ということで、うまく出来たっぽいです。
自前のシミュレーションデータでやってみる
ヘルプの例でうまくいくのは当たり前なので、ちょっと色々小細工したデータでもう一度やってみます。これは以前何度か時系列分析まわりで使ったことのある「周期的スパイク+ベースライン変化を伴う」時系列シミュレーションデータです。
> 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
これもうまくいきました。ではちょっと意地悪してみましょう。
> 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
意外とうまくいってない感がorz という点から察するに、もうちょっと季節性が明確な時系列データでないといけないということなんですかね。。。チャンスがあればちょっと実データにも使ってみようかなと思います。
*1:Johansen's procedureを思い出しますね