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

廿TT

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

ガンマ再生過程に基づくカウントデータの分布

関心のある事象(例えば機械の故障, タクシーの到着など)が繰り返し生起し, それぞれのイベントの生起間隔が独立に同一のガンマ分布に従う場合を考える.

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

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

で表す.

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

f:id:abrahamcow:20170104222627p:plain

上図より,  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)< 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 が正の整数の場合,  F^{(n)}(t) はアーラン分布(アーラン分布 - Wikipedia)になるので,

 F^{(n)}(t) -  F^{(n+1)}(t) は,

 \displaystyle 1-\sum _{k=0}^{n-1}\frac{\lambda x^{k}}{k!}e^{-\lambda x}-(1-\sum _{k=0}^{n}\frac{\lambda x^{k}}{k!}e^{-\lambda x})=\frac{\lambda x^{k}}{k!}e^{-\lambda x}

ポアソン分布になる.

上記の分布はポアソン分布の拡張である.

シミュレーション

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")

f:id:abrahamcow:20170104224340p:plain

付録

ガンマ分布のたたみこみ

\displaystyle f(z) = \int^{\infty}_0 \frac{\beta^k}{\Gamma(k)}x^{k-1}e^{-\beta x}  \frac{\beta^k}{\Gamma(k)}(z-x)^{k-1}e^{-\beta (z-x)} \,dx\\
\displaystyle=\beta^{2k} e^{-\beta z} \frac{1}{2\Gamma(k)} \int^{\infty}_0  x^{k-1}(z-x)^{k-1} \,dx

 x=zu と置くと,  dx=z\,du

\displaystyle f(z) =\beta^{2k} e^{-\beta z}z^{2k-1}  \frac{1}{2\Gamma(k)} \int^{\infty}_0 u^{k-1}(1-u)^{k-1} \, du

あとはベータ関数の性質(ベータ関数 - Wikipedia)より,

\displaystyle f(z) = \frac{\beta^{2k}}{\Gamma(2k)} e^{-\beta z}z^{2k-1}.

ガンマ分布のラプラス変換

ラプラス変換を使うとたたみこみの計算がかんたんになる(場合がある).

 \displaystyle \int^{\infty}_{0} e^{-sx} \frac{\beta^k}{\Gamma(k)}x^{k-1}e^{-\beta x} \,dx\\
= \frac{\beta^k}{\Gamma(k)} \int^{\infty}_{0} x^{k-1} e^{-(s+\beta)x}\,dx

 (\beta+s)x=u と置くと,  du=(\beta+s)\,dx. x=u/(\beta+s)

\displaystyle = \frac{\beta^k}{(\beta+s)^k} \frac{1}{\Gamma(k)} \int^{\infty}_{0} u^{k-1} e^{-u}\,dx
\displaystyle = \frac{\beta^k}{(\beta+s)^k}

ガンマ分布の密度関数を,

\displaystyle f(x) = \frac{1}{\sigma \Gamma(k) x^{k-1}e^{-\frac{x}{\sigma}}} と表す場合もある.

 \beta =\sigma^{-1} と置き換えると, ラプラス変換は,

\displaystyle \frac{1}{(1+\sigma s)^k}.

スバラシク実力がつくと評判のラプラス変換キャンパス・ゼミ

スバラシク実力がつくと評判のラプラス変換キャンパス・ゼミ