廿TT

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

混合ポアソン分布による逐次更新型異常検知をRで

今日の川柳

詳しい説明は『異常検知と変化検知』5章を見てください。

異常検知と変化検知 (機械学習プロフェッショナルシリーズ)

異常検知と変化検知 (機械学習プロフェッショナルシリーズ)

『異常検知と変化検知』では混合正規分布を仮定していますが、ここではポアソン分布でやってみます。

計算メモ

重み付き対数尤度:
「直近の観測は重視するが古い観測は徐々に忘れたい」ため、対数尤度に時刻 t に依存する重み w_t をつける。


 \displaystyle l= \sum_{n=1}^{t} w_t^{(n)} \log \sum_{k=1}^{K} \pi_k p (x^{(n)}|\theta_k)
ここで p(x|\theta)ポアソン分布の確率関数で、 \pi_k は混合比率で、
\displaystyle \sum_{k=1}^{K} \pi_k =1
を満たす。
パラメータの推定は時刻 t までに得られたデータにもとづいておこなう。


和と対数の順序交換:
上記の対数尤度は対数の中に和があるため最大化の計算が面倒。
そこで
\displaystyle \sum_{k=1}^{K} q_k^{(n)} =1
を満たす係数 q_k^{(n)} を導入し、対数尤度を以下のように書き直す。


 \displaystyle l= \sum_{n=1}^{t} w_t^{(n)} \log \sum_{k=1}^{K} q_k^{(n)} \frac{\pi_k p (x^{(n)}|\theta_k)}{q_k^{(n)}}
イェンセンの不等式を使って対数尤度の下限を求める。

 \displaystyle l_{0}= \sum_{n=1}^{t} w_t^{(n)}  \sum_{k=1}^{K} q_k^{(n)} \log \frac{\pi_k p (x^{(n)}|\theta_k)}{q_k^{(n)}}
この下限を最大化することで、対数尤度を最大化する。


q_k^{(n)} の最適化:
まず  \theta_k を既知として q_k^{(n)} を求める。
ラグランジュ乗数 \lambda_n を使って、最適解の条件は


 \displaystyle \frac{\partial}{\partial q_k^{(n)}} = l_{0}-\lambda_n \left(\sum_{k=1}^{K}q_k^{(n)}-1\right)=0
となる。右辺は、

 \displaystyle w_t^{(n)}  \sum_{k=1}^{K} q_k^{(n)} \log \left\{\pi_k p (x^{(n)}|\theta_k)\right\}-(1+\log q_k^{(n)}) -\lambda_n

となり解は、

\displaystyle q_k^{(n)}=\exp\left(-1-\frac{\lambda_n}{w_t^{(n)}}\right)\pi_k p (x^{(n)}|\theta_k)
と得られる。両辺の k に関する和を取り、

\displaystyle 1=\exp\left(-1-\frac{\lambda_n}{w_t^{(n)}}\right)\sum_{k=1}^{K} \pi_k p (x^{(n)}|\theta_k)

であるから、

\displaystyle \exp\left(-1-\frac{\lambda_n}{w_t^{(n)}}\right)=\frac{1}{\sum_{k=1}^{K} \pi_k p (x^{(n)}|\theta_k)}

と求まり、結果、

\displaystyle q_k^{(n)}=\frac{w_t^{(n)}\pi_k p (x^{(n)}|\theta_k)}{\sum_{k=1}^{K} \pi_k p (x^{(n)}|\theta_k)}
となる。


\theta_k の最適化:


\displaystyle  \frac{\partial}{\partial \theta_k}l_{0}= \sum_{n=1}^{t} w_t^{(n)}  q_k^{(n)} \log \{p (x^{(n)}|\theta_k)\} \\
\displaystyle =\sum_{n=1}^{t} w_t^{(n)}  \sum_{k=1}^{K}q_k^{(n)} \left(x^{(n)}/\theta- 1\right)=0
を解く。

\displaystyle \hat \theta_{k}^{(t)}=\frac{1}{\sum_{n=1}^{t} w_t^{(n)} q_k^{(n)}}\sum_{n=1}^{t} w_t^{(n)}q_k^{(n)} x^{(n)}
となる。


\pi_k のスムージング:
事前分布として \pi_k にディリクレ分布を仮定する。ディリクレ分布のパラメータを \alpha_k (k=1,\ldots,K) とし、やはりラグランジュ乗数 \lambda を用いて、


\displaystyle \frac{\partial}{\partial \pi_k}l_{0}=\frac{\partial}{\partial \pi_k}\left[\sum_{n=1}^{t}w_t^{(n)}\sum_{k=1}^{K}q_{k}^{(n)}\log \pi_k + \sum_{k=1}^{K}(\alpha_k-1)\log \pi_k - \lambda \left(\sum_{k=1}^{K}\pi_k-1\right) \right]
を解く。

\displaystyle \sum_{n=1}^{t} w_t^{(n)} q_k^{(n)}+(\alpha_k-1)/\pi_k-\lambda=0
\displaystyle \pi_k=\frac{1}{\lambda}\left( \sum_{n=1}^{t} w_t^{(n)} q_k^{(n)}+(\alpha_k-1)\right)
と得られる。両辺の k に関する和を取り、

\displaystyle 1=\frac{1}{\lambda}\left( \sum_{k=1}^{K}\sum_{n=1}^{t} w_t^{(n)} q_k^{(n)}+\sum_{k=1}^{K}(\alpha_k-1)\right)
\displaystyle \lambda= \sum_{k=1}^{K}\sum_{n=1}^{t} w_t^{(n)} q_k^{(n)}+\sum_{k=1}^{K}(\alpha_k-1)
であるから、

\displaystyle \hat \pi_{k}^{(t)}=\frac{ \sum_{n=1}^{t} w_t^{(n)} q_k^{(n)}+(\alpha_k-1)}{\sum_{k=1}^{K}\sum_{n=1}^{t} w_t^{(n)} q_k^{(n)}+\sum_{k=1}^{K}(\alpha_k-1)}
となる。
\tilde \pi_k^{(t)}=\sum_{n=1}^{t} w_t^{(n)} q_k^{(n)} とおき、\alpha_k=1+\gamma と設定すると、

\displaystyle \hat \pi_{k}^{(t)}=\frac{ \sum_{n=1}^{t} \tilde \pi_k^{(t)}+\gamma}{K\gamma + \sum_{k=1}^{K}\tilde \pi_k^{(t)}}
となる。


重みの選択:


\displaystyle w_t^{(n)}=\beta(1-\beta)^{t-n}
(0<\beta<1)とすると、

\displaystyle \tilde \pi_k^{(t-1)}=\beta\sum_{n=1}^{t-1}(1-\beta)^{t-n-1}q_k^{(n)}
\displaystyle \tilde \pi_k^{t}=\beta\sum_{n=1}^{t}(1-\beta)^{t-n}q_k^{(n)}\\
\displaystyle = \beta q_k^{(t)}+\beta \sum_{n=1}^{t-1}(1-\beta)^{t-n}q_{k}^{(n)}\\
\displaystyle = \beta q_k^{(t)}+(1-\beta)\tilde \pi_k^{(t-1)}
という更新式が成り立つ。
\hat \theta_k については、 \tilde \theta_k = \sum_{n=1}^{t} w_t^{(n)}q_k^{(n)} x^{(n)} とおくと、

\displaystyle \tilde \theta_{k}^{(t)}=\beta q_k^{(t-1)} x^{(t)}+ (1-\beta)\tilde \theta_{k}^{(t-1)}
という更新式が成り立つ。

R による実装例

seqmixpois <-function(y,K,beta,gamma,lambda_est,pi_est){
  Tmax<- length(y)
  lambda_rec <- pi_rec <- matrix(NA,Tmax,K)
  alpha_rec <- numeric(Tmax)
  pi_tilde <- 0
  lambda_tilde <-0
  for (t in 1:Tmax) {
    q_unnorm <- sapply(1:K,function(k)pi_est[k]*dpois(y[t],lambda_est[k]))
    q_norm <- q_unnorm/sum(q_unnorm)
    pi_tilde <- (1-beta)*pi_tilde + beta*q_norm
    lambda_tilde <- (1-beta)*lambda_tilde + beta*q_norm*y[t]
    lambda_est <- lambda_tilde/pi_tilde
    lambda_rec[t,] <- lambda_est
    pi_est <-(pi_tilde+gamma)/(K*gamma+sum(pi_tilde))
    pi_rec[t,] <- pi_est
    alpha_rec[t] <- -log(sum( sapply(1:K,function(k)pi_est[k]*dpois(y[t],lambda_est[k]))))
  }
  list(lambda=lambda_rec,pi=pi_rec,alpha=alpha_rec)
}

適当に乱数を生成して様子を見る。

#パラメータ
pi <- c(0.2,0.3,0.5)
lambda <-c(1,5,10)
Tmax <- 1000
set.seed(1234)
ind_l <-sample.int(3, Tmax, replace = TRUE, prob = pi)
x <- rpois(Tmax,lambda[ind_l])
anomaly <- sample.int(Tmax,3)
x[anomaly] <- rpois(3,30)
pi_ini <- drop(gtools::rdirichlet(1,c(1,1,1)))
lambda_ini <- rexp(3,1/5)
res <-seqmixpois(x,K = 3,beta = 0.1,gamma = 0.05,lambda_est = lambda_ini,pi_est = pi_ini)

matplot(res$pi,type = "l",lty=1,col = c("royalblue","forestgreen","orange2"),ylim=c(0,1),ylab = "")
abline(h=pi,lty=2,lwd=1.5)

matplot(res$lambda,type = "l",lty=1,col = c("royalblue","forestgreen","orange2"),ylab = "")
abline(h=lambda,lty=2,lwd=1.4)

plot(res$alpha,type = "l",ylab = "")
points(anomaly,res$alpha[anomaly],col="red")

f:id:abrahamcow:20180223203422p:plain
\pi_k の推定値。点線が真値。

f:id:abrahamcow:20180223203518p:plain
\theta_k の推定値。真値のまわりをうろうろするようです。

f:id:abrahamcow:20180223203616p:plain
負の対数尤度。大きいほど異常度が高い。赤い丸はわざと別の分布から持ってきた値を入れた点。ただしく異常度が高くなってる。