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

廿TT

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

RStanで離散パラメータを含むモデルの推定(disaster model)

前置き

Stan のマニュアル11章の例題です。
Stan は離散パラメータをサポートしていないので、離散パラメータを含むモデルの推定では周辺化して離散パラメータを消去してやる必要があります。その練習です。

解析対称はイギリスの炭鉱事故の発生件数のデータです。

3. Tutorial — PyMC 2.3.6 documentation から取得しました。

#R のコード
dat <- c(4, 5, 4, 0, 1, 4, 3, 4, 0, 6, 3, 3, 4, 0, 2, 6,
          3, 3, 5, 4, 5, 3, 1, 4, 4, 1, 5, 5, 3, 4, 2, 5,
          2, 2, 3, 4, 2, 1, 3, 2, 2, 1, 1, 1, 1, 3, 0, 0,
          1, 0, 1, 1, 0, 0, 3, 1, 0, 3, 2, 2, 0, 1, 1, 1,
          0, 1, 0, 1, 0, 0, 0, 2, 1, 0, 0, 0, 1, 1, 0, 2,
          3, 3, 1, 1, 2, 1, 1, 1, 1, 2, 4, 2, 0, 0, 1, 4,
          0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 1)
plot(dat,type="h",ylab="frequency")

f:id:abrahamcow:20151029003418p:plain

40年ごろを境に事故の発生件数が減少していることが伺えるので、その変化を抽出できるようなモデルを考えます。

事故の発生件数はポアソン分布に従うと仮定し、ある時期を境にポアソン分布の平均が変化したと考えると良さそうです。

そこで t 年目の事故の発生件数を D_t とし、D_t はパラメータ \lambda _tポアソン分布に従うというモデルを立てました。

ただし \lambda _tt < s のとき e、そうでなければ l の値を取ります。el はそれぞれ early、late の頭文字です。

無情報事前分布として、パラメータ el には区間 [0, 100] の連続一様分布を仮定し、s に離散一様分布を仮定することにします。

確率分布を列挙すると以下のようになります。

\displaystyle e  \sim {\rm continuous~uniform}(0,100)
\displaystyle l  \sim {\rm continuous~uniform}(0,100)
\displaystyle s \sim {\rm discrete~uniform}(1,T)
\displaystyle D_t \sim {\rm Poisson}(\lambda_t)

同時確率を書き下すと

\displaystyle p(e) p(l) p(s) p(D|s,e,l)

となります。

周辺化

 p(e)p(l) は良いとして、s は離散パラメータなので周辺化します。

p(s) = 1/T, ~(s=1,\ldots,T) なので、

 p(D|e, l) \\
= \sum_{s=1}^{T} p(s, D|e, l) \\
= \sum_{s=1}^{T} p(s) p(D|e, l) \\
= \sum_{s=1}^{T}\left( (1/T) \prod_{t=1}^{T} {\rm Poisson}(D_t|\lambda_t) \right)

対数をとり、

 \log p(D|e, l) = \\
\log \sum_{s=1}^{T}\left( (1/T) \prod_{t=1}^{T} {\rm Poisson}(D_t|\lambda_t) \right) \\
=\log \sum_{s=1}^{T} \exp \left( \log(1/T)+\sum_{t=1}^{T} \log {\rm Poisson}(D_t|\lambda_t) \right)

コーディング

#Stan のコード
data {
      int<lower=1> T;
      int<lower=0> D[T];
}
transformed data {
      real log_unif;
      log_unif <- -log(T);
}
parameters {
      real<lower=0> e;
      real<lower=0> l;
}
transformed parameters {
      vector[T] lp;
      lp <- rep_vector(log_unif, T);
      for (s in 1:T){
      	 for (t in 1:T){
      	 	lp[s] <- lp[s] + poisson_log(D[t], if_else(t < s, e, l));
      	 }
      }
} model {
      e ~ uniform(0,100);
      l ~ uniform(0,100);
      increment_log_prob(log_sum_exp(lp));
}
 generated quantities {
      int<lower=1,upper=T> s;
      s <- categorical_rng(softmax(lp));
}

 \log \sum \exp() は Stan では log_sum_exp という関数で表現されます。

e ~ uniform(0,100)、 l ~ uniform(0,100) は \log p(e)\log p(l)対数事後密度に足し算するという意味を持っています。

パラメータ s の分布も気になるので generated quantities ブロックで乱数を生成しています。

softmax(lp) は

\displaystyle \frac{\exp(\log(1/T))}{\sum\exp(\log(1/T))}

の意味です。

結果

パラメータ s の分布は以下のようになりました。

f:id:abrahamcow:20151029004039p:plain

EAP 推定値 \hat \lambda _t をデータに重ねてプロットしたのが下の図です。

f:id:abrahamcow:20151029004904p:plain

サンプリングの結果を重ね書きしてパラメータの分布を見たのが下の図です。

f:id:abrahamcow:20151029004748p:plain

# R のコード
library(rstan)
dat4stan <-list(T=length(dat),D=dat)
cpm <- stan_model("changepoint.stan")
fit <- sampling(cpm, data=dat4stan)
#traceplot(fit,pars=c("e","l"))
ex <- extract(fit)
plot(table(ex$s),ylab="")
ests <-mean(ex$s)
este <-mean(ex$e)
estl <-mean(ex$l)
lambda_t <- function(t){
  ifelse(t>ests,estl,este)
}
plot(dat,type="h",ylab="frequency")
curve(lambda_t(x),add=TRUE,col="royalblue",lty=2,lwd=2)

plot(dat,type="n",ylab="")
for(i in 1:4000){
  ests <-ex$s[i]
  este <-ex$e[i]
  estl <-ex$l[i]
  curve(lambda_t(x),add=TRUE,col="cornflowerblue")
}

おまけ

パラメータの推定をうまくいってることを確かめるには、以下のように乱数でデータセットを生成してやるとよいでしょう。

e <- 7
l <- 2
s <-  20
T <- 50
lambda <-c(rep(e,s),rep(l,T-s))
set.seed(1)
D <-rpois(T,lambda)
plot(D,type="h")
dat4stan <-list(T=T,D=D)

参考文献と関連エントリ

Stan - Documentation

岩波データサイエンス Vol.1

岩波データサイエンス Vol.1

基礎からのベイズ統計学: ハミルトニアンモンテカルロ法による実践的入門

基礎からのベイズ統計学: ハミルトニアンモンテカルロ法による実践的入門

abrahamcow.hatenablog.com