ジョンとヨーコのイマジン日記

library(estatapi) library(tidyverse) myappId <- scan("appId.txt", what = character()) #ここには自分のアプリケーションIDを入れる #### dat1 <-estat_getStatsData(appId = myappId, statsDataId = "0000010209") unique(dat1$`Ｉ　健康・医療`) nyuin <-dplyr::filter(dat1,`Ｉ　健康・医療`=="#I04103_精神科病院年間新入院患者数（人口10万人当たり）" , 地域=="全国") %>% mutate(year=as.integer(substr(調査年,1,4))) theme_set(theme_bw(12,"Osaka")+ theme(axis.text = element_text(colour="black"))) p_nyuin <-ggplot(nyuin,aes(x=year,y=value))+ geom_line()+ labs(y="人口10万人当たり", x="年", title="精神科病院年間新入院患者数（全国）\n社会・人口統計体系都道府県データ / 社会生活統計指標より") print(p_nyuin) ggsave("~/Desktop/p_nyuin.png",p_nyuin, width = 7, height = 7) #### zaiin <-dplyr::filter(dat1,`Ｉ　健康・医療`=="#I10205_精神科病院平均在院日数" , 地域=="全国") %>% mutate(year=as.integer(substr(調査年,1,4))) p_zaiin <-ggplot(zaiin,aes(x=year,y=value))+ geom_line()+ labs(y="人口10万人当たり", x="年", title="精神科病院平均在院日数（全国）\n社会・人口統計体系都道府県データ / 社会生活統計指標より") print(p_zaiin) ggsave("~/Desktop/p_zaiin.png",p_zaiin, width = 7, height = 7) #### byosho <-dplyr::filter(dat1,`Ｉ　健康・医療`=="#I0910205_精神病床数（人口10万人当たり）", 地域=="全国") %>% mutate(year=as.integer(substr(調査年,1,4))) p_byosho <-ggplot(byosho,aes(x=year,y=value))+ geom_line()+ labs(y="人口10万人当たり",x="年", title="精神病床数（全国）\n社会・人口統計体系都道府県データ / 社会生活統計指標より") print(p_byosho) ggsave("~/Desktop/p_byosho.png",p_byosho, width = 7, height = 7)

再生過程とは関心のある事象（例えば機械の故障，タクシーの到着など）が繰り返し生起し，それぞれのイベントの生起間隔が独立に同一の分布 $F(t)$ に従う確率過程である．

今回は $F(t)$ がガンマ分布の場合を考える．

イベントの生起間隔を確率変数 $t_i$ で表す. またイベントの発生時刻は，

$T_n=\sum_{i=1}^{n}t_i$

で表す.

いま, 区間 $(0,t)$ で起こったイベントの数を $N(t)$ で表し, この $N(t)$ の分布が知りたい.

上図より, $N(t) \lt n$ と $T_n \ge t$ は同値なので,

$P (N(t)\lt n) = P (T_n \ge t)$ .

$P (T_n \ge t)$ を求めるには, ガンマ分布の n 重たたみこみ $F^{(n)}(t)$ を求めればよい. ガンマ分布の再生性より，

$\displaystyle F^{(n)}(t) = \frac{1}{\Gamma(nk)} \int ^{\beta t} _{0} u^{n k-1} e^{-u} \, du$ .

よって，求めたい確率は，

$P(N(t)= n)= P(N(t) \lt n+1) - P(N(t) \lt n)\\ =1-F^{(n+1)}(t) - (1- F^{(n)}(t))\\ = F^{(n)}(t) - F^{(n+1)}(t)$

これを R で計算するには，

pgamma(t,shape*k,rate)-pgamma(t,shape*k+shape,rate)

とすればよい.

形状パラメータ k が正の整数の場合はポアソン分布になる.

シミュレーションしてみる：

iter =10000
Ns <- numeric(iter)
Tim <- 10
for(i in 1:iter){
  ti <-cumsum(rgamma(100,1.5,2))
  Ns[i] = sum(ti<Tim)
}
Gpmf <- function(k,shape,rate,t){
  pgamma(t,shape*k,rate)-pgamma(t,shape*k+shape,rate)
}
minN<-min(Ns)
maxN<-max(Ns)
pred <- Gpmf(minN:maxN,1.5,2,Tim)*iter 
plot(table(Ns),xlab="",ylab="")
points(minN:maxN,pred,type="b")

参考文献はこちら：

Winkelmann, R. (1995). Duration dependence and dispersion in count-data models. Journal of business & economic statistics, 13(4), 467-474. http://www.econ.uzh.ch/dam/jcr:bb267be8-483a-45d1-90d3-338016548f60/jbes.pdf

ジョンとヨーコのイマジン日記

キョウとアンナのラヴラブダイアイリー改め、ジョンとヨーコのイマジン日記です。

データに見る精神疾患の軽症化とびまん化（あるいはe-tstat APIのデモ）

ガンマ分布の再生過程における再生回数の分布