廿TT

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

状態空間モデル(ローカルレベルモデル)で変化点の検出

ローカルレベルモデル

ローカルレベルモデルは以下の組で表されるモデルです.

観測モデル:
\displaystyle y_t = \alpha_t + \varepsilon_t
\varepsilon_t \sim {\rm Normal}(0,\sigma^2_\varepsilon)

システムモデル:
\displaystyle \alpha_t = \alpha_{t-1} + \eta_t
\eta_t \sim {\rm Normal}(0,\sigma^2_\eta)

y_t は観測値, \alpha_t は未観測の「状態」です.

\varepsilon_t は観測ノイズ, \eta_t はシステムノイズと呼ばれます.

\displaystyle \alpha_t = \alpha_{t-1} + \eta_t
という式は「状態は変化しているが, 一歩前の値と近い値を取るだろう」という仮説を反映しています.

\displaystyle \alpha_t \approx \alpha_{t-1}

裏を返せば, \alpha_t が大きく変化したところ, すなわち \eta_t が大きいところは, 状態が変化したポイントです.

dlm

R の dlm パッケージで  \hat \eta_t を求めてみます.

library(dlm)
N<-length(Nile)
build.1 <- function(theta){
  dlmModPoly(order=1, dV=exp(theta[1]), dW=exp(theta[2]))
}
fit.1 <- dlmMLE(Nile, parm=c(1, 1), build.1)
mod.Nile <- build.1(fit.1$par)
NileFilt <- dlmFilter(Nile, mod.Nile)
NileSmooth <- dlmSmooth(NileFilt)
etahat <-diff(NileSmooth$s[-1])
x<-(start(Nile)[1]:end(Nile)[1])
cp <-which.max(abs(etahat))

変化点は28番目にありました.

> cp
[1] 28

図示.

f:id:abrahamcow:20171229212900p:plain

KFAS

KFASパッケージでも同じことができます.

KFASで求めた  \hat \eta_t とdlmで求めた  \hat \eta_t はだいたい近い値になりました.

library(KFAS)
mod2 <- SSModel(Nile~SSMtrend(degree = 1, Q = NA),H=NA)
fit2 <- fitSSM(mod2,numeric(2))
kfs2 <-KFS(fit2$model,smoothing = c("state","mean","disturbance"))
ran1 <-range(c(kfs2$etahat[-N],etahat))
plot(kfs2$etahat[-N],etahat,xlim=ran1,ylim=ran1)
abline(0,1,lty=2)

f:id:abrahamcow:20171229213115p:plain

結果は一緒です.

> which.max(abs(kfs2$etahat))
[1] 28

KFASでは関数rstandardを使うと標準化したシステムノイズを返してくれるので, 検定の枠組みに持っていくこともたぶんできます.

でも実用的には  \hat \eta_t閾値を超えた点を変化点とみなすことにしたらいいかなと思います.

stdresid <-rstandard(kfs2,"state")
hist(stdresid[-N],freq = FALSE)
curve(dnorm(x),add=TRUE,lwd=2)

f:id:abrahamcow:20171229213506p:plain

changepoint

変化点を検出するパッケージにchangepointというのがあります.

これを使っても結果は一緒になりました.

> library(changepoint)
> ans <-cpt.mean(Nile)
> summary(ans)
Created Using changepoint version 2.2.2 
Changepoint type      : Change in mean 
Method of analysis    : AMOC 
Test Statistic  : Normal 
Type of penalty       : MBIC with value, 13.81551 
Minimum Segment Length : 1 
Maximum no. of cpts   : 1 
Changepoint Locations : 28

cusum管理図

cusum管理図というのもあります.

同じ結果です.

> cusum1 <-cumsum(Nile-mean(Nile))
> which.max(cusum1)
[1] 28