読者です 読者をやめる 読者になる 読者になる

廿TT

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

状態空間非定常ポアソン(NHPP using Stan)

R 確率過程 Stan

ポアソン過程は再発事象のモデルとしてよく使われる。

ポアソン過程ではイベントが観測された時刻を T_t (t=1,\ldots,N) とすると、イベントの生起間隔  \Delta T_t = T_t - T_{t-1} は独立にパラメータ λ の指数分布に従う。

ポアソン過程の拡張としてパラメータλ が時間に依存して変化する非定常ポアソン過程というのが考えられる。

今回は以下のように一回差分のトレンドで時間変化を表現する。

 \mu_t \sim {\rm Normal}(\mu_{t-1},\sigma^2),~t=2,\ldots,N
 \Delta T_t \sim {\rm Exponential}(\exp(\mu_t)) ~ t=2,\ldots,N

\mu_1\sigma には幅の広い一様分布を仮定する。

Stan コードはこうなりました。

//NHPP.stan
data {
  int<lower=1> n;
  real<lower=0> D[n];
}
parameters {
  vector[n] loglambda;
  real<lower=0>sigma;
}
model {
  loglambda[2:n] ~ normal(loglambda[1:(n-1)],sigma);
  D ~ exponential(exp(loglambda));
}

R の boot パッケージには炭鉱災害の発生を記録した coal データというのが入っている。

data("coal",package = "boot")
p1 <-ggplot(coal)+
  geom_step(aes(x=date,y=1:length(date)))+
  ylab("Cumulative recurrences")
print(p1)

f:id:abrahamcow:20161126085429p:plain

これに上記の NHPP(Non-Homogeneous Poisson Process; 非定常ポアソン過程)を当てはめる。

rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())
NHPPmodel <-stan_model("NHPP.stan")
dat <-list(D=diff(ti),n=n-1)
NHPPfit <-sampling(NHPPmodel,dat,iter=10000)
lambdahat <- exp(get_posterior_mean(NHPPfit,pars="loglambda")[,"mean-all chains"])
exNHPP <-rstan::extract(NHPPfit)
lower <-exp(apply(exNHPP$loglambda, 2, quantile,0.2))
upper <-exp(apply(exNHPP$loglambda, 2, quantile,0.8))
p1+geom_line(aes(x=cumsum(c(date[1],1/lambdahat)),y=1:length(date)),colour="royalblue",size=1)+
  geom_line(aes(x=cumsum(c(date[1],1/lower)),y=1:length(date)),colour="royalblue",size=1,linetype=2)+
  geom_line(aes(x=cumsum(c(date[1],1/upper)),y=1:length(date)),colour="royalblue",size=1,linetype=2)

f:id:abrahamcow:20161126085730p:plain

点線は60%のベイズ信頼区間です。

これがなんの役に立つのかは後で考える。

2回差分のトレンドでもやってみました。

data {
  int<lower=1> n;
  real<lower=0> D[n];
}
parameters {
  vector[n] loglambda;
  real<lower=0>sigma;
}
model {
  loglambda[3:n] ~ normal(2*loglambda[2:(n-1)]-loglambda[1:(n-2)],sigma);
  D ~ exponential(exp(loglambda));
  sigma ~ exponential(0.1);
}

sigma が収束しにくかったので指数分布を入れて縛ってます。

f:id:abrahamcow:20161126090235p:plain

そのせいかベイズ信頼区間の幅はちょっと狭くなった。

グラフを見た限りでは、インテンシティ(loglambda)の経時的な変化はどちらでも大差ない感じ。

StanとRでベイズ統計モデリング (Wonderful R)

StanとRでベイズ統計モデリング (Wonderful R)

abrahamcow.hatenablog.com