R 時系列データの異常部位検出（近傍法） - 統計学と疫学と時々、助教生活

今回も先ほど公開した記事と同様に、『入門機械学習による異常検知-Rによる実践ガイド-』（コロナ社、井出剛著、2015）の中にある「7.1 近傍法による異常検知」ということで、この内容をなぞり、最後に時系列データ解析の際に重要になる窓幅（どれだけの地点のデータをまとめるか。年ごとのデータで窓幅が3といえば、3年分の値をまとめてベクトルにするイメージを持てば良い。）を変化させると、どのようになるかも図示する。この本で取り上げているのは、Keoghらの心電図のデータであり、今回もこれをhttp://www.cs.ucr.edu/~eamonn/discords/qtdbsel102.txtからダウンロードする（2017/4/28現在ダウンロード可）。

時系列データの異常検知とは（P.196）

この場合は、窓と呼ばれる単位でいくつかのデータを一つのベクトル（「部分時系列」と本書では呼んでいる）にまとめ上げる。この場合、他の部分時系列とは異なる振る舞いをする部分を同定することを指す。このことを本書では異常部位検出と呼んでいるが、外れ値検出問題を時系列に当てはめたと考えることができるため、異常度には近傍法を用いることができる。

データの読み込みと図示

X <- read.table("qtdbsel102.txt")
Xtr <- X[1:3000,2] #訓練データ
X <- X[3001:6000,2] #検証データ
plot(Xtr, type="l", ylim=c(4.0,6.0)) #訓練データの図示
plot(X, type="l", ylim=c(4.0,6.0)) #検証データの図示

訓練用データ

f:id:ryosukefujii0320:20170428220610p:plain

検証用データ

f:id:ryosukefujii0320:20170428220640p:plain

こんな感じでデータを2分割して、訓練用データの各要素と検証用データの各要素の距離を計算する。

窓幅を100をデフォルトにして、最近傍までの距離を異常度にして行う。

X <- read.table("qtdbsel102.txt")
w <- 100; nk <- 1 #窓の数と近傍数
Xtr <- X[1:3000,2]; Dtr <- embed(Xtr, w) #訓練データ
X <- X[3001:6000,2]; D <- embed(X,w) #検証データ
d <- knnx.dist(Dtr, D, k=nk); a <- d[,1] #最近傍までの距離
plot(a, ylab="annomaly score", type="l") #結果の図示

f:id:ryosukefujii0320:20170428221433p:plain

窓幅を10, 50, 200, 1000と変化させて図示する。

X <- read.table("qtdbsel102.txt")
window <- c(10, 50, 200, 500, 1000) 
for(w in window){
  Xtr <- X[1:3000,2]; Dtr <- embed(Xtr, w)
  Xval <- X[3001:6000,2]; Dval <- embed(Xval,w) 
  d <- knnx.dist(Dtr, Dval, k=nk); a <- d[,1] 
  plot(a, main=paste("width of window =", w), ylab="annomaly score", type="l") 
}