廿TT

譬如水怙牛過窓櫺 頭角四蹄都過了 因甚麼尾巴過不得

カルバック・ライブラー密度比推定法をRで(異常検知と変化検知)

『異常検知と変化検知』の12章で説明されているカルバック・ライブラー密度比推定法をRで試してみました.

正しく実装できてる自信がないので, コメントいただけると嬉しいです.

アルゴリズム

アルゴリズムは以下のような感じだと思っています.

正常標本とみなす訓練データを \boldsymbol{x}=(x_1,\ldots,x_N)^T , 異常標本が含まれるかもしれないテストデータを \boldsymbol{y}=(y_1,\ldots,y_M)^T とする.

バンド幅 h, 勾配法のステップ幅 \eta, \boldsymbol{\theta}=(\theta_1,\ldots, \theta_n)^T の初期値を与える.

\displaystyle \psi_n(u) = \exp\left(-\frac{\sqrt{(u-x_n)^2}}{2h^2}\right)
とする(他のカーネルを用いてもよい).

\boldsymbol{\psi}(u) = (\psi_1(u),\dots,\psi_N(u))^T とする.

\displaystyle \nabla J(\boldsymbol{\theta}) = \frac{1}{M}\sum_{m=1}^{M}\boldsymbol{\psi}(y_m) - \frac{1}{N}\sum_{n=1}^{N}\frac{\boldsymbol{\psi}(x_n)}{\boldsymbol{\theta}^T \boldsymbol{\psi}(x_n)}
とする.

\displaystyle \boldsymbol{\theta} \leftarrow \boldsymbol{\theta} -\eta \nabla J(\boldsymbol{\theta})
を更新式として \boldsymbol{\theta} が収束するまで繰り返す.

\displaystyle r_\boldsymbol{\theta}(u) =\sum_{n=1}^{N}\theta_n\exp\left(-\frac{\sqrt{(u-x_n)^2}}{2h^2}\right)
を密度比の推定量とする.

\displaystyle -\log r_{\boldsymbol{\theta}}(y_m)
を各標本の異常度とする.

なぜこうなるのかについては本を見てください.

R による実装例

KLdre <-function(train,test,ini,eta=0.01,tol=1e-4,maxit=10000){
  RBF_kernel <-function(x,y,bw){
    exp(-(sqrt(((x-y))^2))/(2*bw^2))
  }
  N <- length(train)
  h=bw.nrd(train)
  psi_train <- t(sapply(train,function(x)RBF_kernel(x,train,h)))
  psi_test <- t(sapply(test,function(x)RBF_kernel(x,train,h)))
  psi_test_mean <- apply(psi_test,2,mean)
  theta <- ini
  for(i in 1:maxit){
    den <- drop(t(theta) %*% psi_train)
    DeltaJ <- psi_test_mean - apply(sweep(psi_train,2,den,"/"),2,mean)
    theta_new <- theta - eta*DeltaJ
    theta_new<-ifelse(theta_new<0,0, theta_new)
    if(all(abs(theta_new-theta)<tol)){
      break
    }
    theta <- theta_new
  }
  rhat <- apply(sweep(psi_test,2,theta_new,"*"),1,sum)
  score <- -log(rhat)
  list(result=data.frame(number=1:length(test),test,rhat,score),theta=theta_new,iter=i)
}

\theta は(たぶん)0より大きくないといけないので, 0を下回った場合は無理やり0と置いてます.

sapply とか apply とか sweep のところ, あってる自信がまったくない. 雰囲気で書いてます.

バンド幅はクロスバリデーションで決めることが推奨されているのですが, ここでは bw.nrd という関数で決めました. bw.nrd の方法について, 詳しくは Scott, D. W. (1992) Multivariate Density Estimation: Theory, Practice, and Visualization. Wiley. に書いてあるそうです.

標準正規分布から100個乱数を生成して, これを訓練データとします.

また, 独立に標準正規分布から50個乱数を生成して, 10番目を5, 45番目を-5という異常値に入れ替えます.

train <- rnorm(100)
test <- rnorm(50)
test[c(10,45)] <- c(-5,5)

このデータで密度比を推定します.

パラメータが多くてちゃんと推定できるか不安なので, 初期値を乱数で振って2回試してみます.

dr1 <- KLdre(train,test,ini = runif(100))
dr2 <- KLdre(train,test,ini = runif(100))
plot(dr1$theta,dr2$theta)
abline(0,1,lty=2)

\boldsymbol{\theta} に関しては全然安定してません.

plot(dr1$theta,dr2$theta)
abline(0,1,lty=2)

f:id:abrahamcow:20171228175325p:plain

でも密度比に関してはだいたい同じような値が求まったので, まあいいかなと思いました.

plot(dr1$result$rhat,dr2$result$rhat)
abline(0,1,lty=2)

f:id:abrahamcow:20171228175423p:plain

異常度のプロットです.

library(ggplot2)
ggplot(dr1$result,aes(x=number,y=score))+
  geom_label(aes(label=number))

f:id:abrahamcow:20171228175536p:plain

10番目と45番目が正しく異常度が高いです.

密度比のプロットです.

ggplot(dr1$result,aes(x=number,y=rhat))+
  geom_hline(yintercept = 0,linetype=2)+
  geom_label(aes(label=number))

f:id:abrahamcow:20171228175732p:plain

10番目と45番目が0に近いです.

元のデータのプロットです.

ggplot(dr1$result,aes(x=number,y=test))+
  geom_label(aes(label=number))

f:id:abrahamcow:20171228175917p:plain

単変量だったら目視でも外れ値がわかりますね.

感想

  • set.seed するのを忘れてた
  • 収束が遅い. R 向きではないかも
  • 異常度の閾値ってどう決めたらいいのかな