廿TT

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

変分ベイズを使って変化点検知をしてみる(ポアソン過程版)

今日の川柳

記法

イベントの発生時刻を t_i として、\Delta t_i = t_i-t_{i-1} とする。

ぜんぶで N+1 回のイベントが観測されたとする。

また最初のイベントの生起をポアソン過程がはじまった時刻として t_0 とする。

モデル

ポアソン過程ではイベントが生起してから次のイベントまでの待ち時間は指数分布に従う。

 \Delta t_i \sim \mathrm{Exponential}(\lambda_i)

ただし、
i<\tau のとき \lambda_i = \lambda^{(1)}
i\ge\tau のとき \lambda_i = \lambda^{(2)}
とする。

\tau がどこにあるのかはわからないので \tau の事前分布は離散一様分布とする。
 \lambda^{(j)} ( j=1,2) の事前分布はガンマ分布とする。

変分推論

machine-learning.hatenablog.com

アルゴリズムの導出は上の記事とほぼ一緒です。

実装

R です。

VIcp <- function(y,a=1,b=1,iter=10){
  Delta_t <- diff(y)
  N <- length(Delta_t)
  lambda <- rgamma(2,1,1)
  loglambda <- log(lambda)
  for(i in 1:iter){
    lp1 <- loglambda[1]-lambda[1]*Delta_t 
    lp2 <- loglambda[2]-lambda[2]*Delta_t
    lp <- sapply(1:(N-1), function(i)sum(lp1[1:i])+sum(lp2[(i+1):N]))
    lp <- c(lp,sum(lp1))
    pi <- exp(lp)/sum(exp(lp))
    Es1 <- c(rev(cumsum(rev(pi[-1]))),0)
    Es2 <- cumsum(pi)
    ahat1 <- sum(Es1)+a
    ahat2 <- sum(Es2)+a
    bhat1 <- sum(Es1*Delta_t)+b
    bhat2 <- sum(Es2*Delta_t)+b
    lambda <- c(ahat1/bhat1,ahat2/bhat2)
    loglambda <- c(digamma(ahat1)-log(bhat1), digamma(ahat2)-log(bhat2))
  }
  return(list(lambda=lambda,pi=pi))
}

例題

炭鉱災害のデータ(coal)を使って試してみます。

炭鉱災害の累積発生件数です。

f:id:abrahamcow:20190607214230p:plain

変化点 \tau のありそうな位置です。

f:id:abrahamcow:20190607214259p:plain

一年間の平均事故発生件数です。(変化前と変化後)

> print(out$lambda)
[1] 3.1076971 0.9390008

コードです。

library(ggplot2)
data("coal",package = "boot")
out <- VIcp(coal$date)
print(out$lambda)
theme_set(theme_bw(24))
qplot(x=coal$date,y=1:nrow(coal),geom="step")+
  xlab("")+ylab("")
qplot(x=1:(nrow(coal)-1),y=out$pi,geom="col",width=1)+
  xlab("")+ylab("")