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

廿TT

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

ポアソン過程の定義. ポアソン分布と指数分布の関係.

定義

事象が発生するまでの待機時間 X_1,X_2,\ldots, は独立同分布で各確率変数はパラメータ λ の指数分布に従うとする.

また  n \ge 1 に対して  S_n = X_1 + \cdots + X_n かつ  S_0 =0 とする.

このとき  N(t) = \max\{n: S_n \le t \} をパラメータ λポアソン過程と呼ぶ.

分布

さて, N(t) が指数過程でなくポアソン過程と呼ばれる理由を納得するために, N(t) の分布を求める.

場合 1:n = 0 のとき

N(t) = 0 とは, 最初の事象の生起までの待ち時間が, t よりも長いということなので,

\displaystyle P(N(t) = 0) = P(X_1 >t) = e^{\lambda t}

場合 2:n自然数のとき

f:id:abrahamcow:20150620133848p:plain

まず  N(t) \ge nS_n \le t が同値であることに注意する.

\displaystyle P(N(t) = n) = P(N(t) \ge n) - P(N(t) \ge n + 1) \\
\displaystyle  = P(S_n \le t) - P(S_{n+1} \le t)

 S_n = X_1 + \cdots + X_n はパラメータ λ, n のガンマ分布に従う(ときわ台学/統計学/ガンマ分布と指数分布 )ので,

\displaystyle P(N(t) = n) = \int^{t}_0 \frac{\lambda ^n x^{n-1}}{(n-1)!}e^{-\lambda x} \, dx - \int^{t}_0 \frac{\lambda ^{n+1} x^{n}}{n!}e^{-\lambda x} \, dx \\
\displaystyle = \int^{t}_0 \frac{\lambda ^n x^{n-1}}{(n-1)!}e^{-\lambda x}  - \frac{\lambda ^{n+1} x^{n}}{n!}e^{-\lambda x} \, dx \\
\displaystyle = \left[ \frac{(\lambda x)^n}{n!}e^{-\lambda x} \right]^{t}_{0} \\
\displaystyle = \frac{(\lambda t)^n}{n!}e^{-\lambda t}

これはパラメータ λtポアソン分布の確率関数である.

R によるシミュレーション

待機時間が指数分布に従うとき 1 単位時間あたりの事象の発生数はポアソン分布に従う.

f:id:abrahamcow:20150726103615p:plain

N <- 10000
rate1 <- 3

res <- numeric(N)
for(i in 1:N){
  t1<-cumsum(rexp(100, rate1))
  res[i] <-sum(t1 <= 1)
}

tab <-table(res)

barplot(rbind(table(res), dpois(0:max(res), rate1)*N),
        col = c("orange", "cornflowerblue"),
        legend.text = c("Observed", "Theoritical"),
        args.legend = list(x = "topright"),
        beside = TRUE
)

待機時間が他の分布, 例えばワイブル分布に従うとき, 事象の発生数をポアソン分布で記述しようとしてもうまくいかない.

f:id:abrahamcow:20150726104331p:plain

shp <- 3
scl <- 1/5

res <- numeric(N)
for(i in 1:N){
  t1<-cumsum(rweibull(100, shp, scl))
  res[i] <-sum(t1 <= 1)
}

tab <-table(res)
lambda_hat <-mean(res)

barplot(rbind(tab, dpois(min(res):max(res), lambda_hat)*N),
        col = c("orange", "cornflowerblue"),
        legend.text = c("Observed", "Theoritical"),
        args.legend = list(x = "topright"),
        beside = TRUE
)

関連エントリ

abrahamcow.hatenablog.com