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

廿TT

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

RStan でポアソン過程の停止時刻を推定する

R Stan 確率過程

問題設定

n = 20 人の患者さんが繰り返し訪問する病院の窓口を T = 100 時間観察した。

患者 i はそれぞれある時期 t_i に達すると病院への訪問をやめる(理由は不明:病気が治ったか、遠くに引っ越したか、通院に飽きたか、死んだか)。

この訪問がレート λ (全患者共通のパラメータ)のポアソン過程であるとして、各患者さんが訪問をやめた時間 t_i を求めたい。

ただし、各患者さんの訪問回数 x_i は記録されているが、訪問間隔は記録されていない。

シミュレーション

まず乱数を使って擬似的なデータセットを生成する。

ポアソン過程なので到着間隔は指数分布に従う。ここでは λ = 1 とした。

停止時刻 t_i は一様乱数で決めた。

#Rのコード
n <- 20
maxT <- 100
t1list <- vector("list",n)
for(i in 1:20){
  tmp <- cumsum(rexp(100)) #多めに出しとく
  t1list[[i]] <- tmp[tmp<=maxT]
}
maxt1 <- sapply(t1list,max)
t2 <- runif(n,0,maxt1)
x <- mapply(function(x,y){sum(x<y)},t1list,t2)
dat4stan<- list(x=x,n=n, maxT = maxT)

生成されたデータセットはこんな感じ。

> dat4stan
$x
 [1] 11 94 76 68  9  1 86 27 27 53 92 32  3 85 12 47 73 23 69 25

$n
[1] 20

$maxT
[1] 100

推定

λt_i を推定するための Stan コードはこうした。

#"poisson.stan"
data { 
int<lower=0> n;
real<lower=0> maxT;
real<lower=0> x[n];
}
parameters {
	real<lower=0> lambda;
	real<lower=0,upper=maxT> T[n];
}
model {
	for(i in 1:n){
	increment_log_prob(x[i]*log(lambda*T[i])-lambda*T[i]);
	}
}

ポアソン過程の対数尤度 x_i\log(\lambda t_i )-\lambda t_i を increment_log_prob の中に書いた。

#R のコード
library(rstan)
pois <- stan_model(file="poisson.stan")
fit <- sampling(pois, data=dat4stan)

推定結果

MCMC ステップのトレースプロット。

traceplot(fit,pars="lambda")
ex <- extract(fit)

f:id:abrahamcow:20150912142400p:plain

λ は収束してくれたっぽい。

λヒストグラムもかいてみる。

hist(ex$lambda)
estlambda <- mean(ex$lambda)
abline(v=estlambda,col="tomato",lwd=2)
abline(v=1,col="royalblue",lwd=2)

f:id:abrahamcow:20150912142832p:plain

赤い線が推定値(EAP)、青い線が真値。けっこうバイアスがありそう。

> estlambda
[1] 0.9839703

次に t_i の真値と推定値を比較する。

#棒グラフ
estT <-apply(ex$T,2,mean)
barplot(rbind(t2,estT),
        xlab="t",
        col = c("orange","cornflowerblue"),
        legend.text = c("true", "estimates"),
        args.legend = list(x = "topright"),
        ylim=c(0,100),
        beside = TRUE)

f:id:abrahamcow:20150912143536p:plain

plot(t2,estT,xlim=c(0,100),ylim=c(0,100))
abline(0,1,lwd=1.5,lty=2)

f:id:abrahamcow:20150912143858p:plain

縦軸を推定値、横軸を真値にしたプロット。

もし完全に一致していたら斜めの線上にぴたっと並ぶ。

まあよく求まっているんじゃないかと思う。

CI <- apply(ex$T,2,function(x)quantile(x,c(0.025, 0.975)))
library(ggplot2)
library(pipeR)
library(tidyr)
ggplot(data.frame(id=1:n,t2=t2),aes(x=id,y=t2)) +
  geom_point(size=3,colour="royalblue") +
  geom_point(size=3,aes(x=id,y=estT),colour="tomato") +
  geom_errorbar(aes(x=1:n,ymin=CI[2,],ymax=CI[1,]),colour="tomato") +
  theme_grey(20)

f:id:abrahamcow:20150912144511p:plain

95%信用区間のエラーバー。青い点が真値、赤が推定値。