私の今後は？潜伏期・再発 - ryamadaの遺伝学・遺伝統計学メモでガンマ分布を扱っていた。ガンマ分布について調べていると指数分布にも関係があるようだったのでメモ。
1. 特徴
指数分布
指数分布はあるランダムなイベントが起こるまでの時間の間隔に対する分布(生起期間の確率)。$P(X > s+t| X > s) = P(X > t)$を満たし、無記憶性がある。つまり、事象が起こったかどうかは次の区間の確率に影響しない。 銀行窓口に客が来る間隔，機械が故障する間隔、死亡するまでの時間*1、電話の待ち時間などがこれに従うと言われている。一般教養の先生は「A君に彼女ができる確率が指数分布だったら」みたいな話をしていた気がする。
ガンマ分布
指数分布は一度イベントが一度起こるまでの間隔を表すのに対し、ガンマ分布は単位時間あたりの回数が$\dfrac{1}{\beta}$のイベントが$\alpha$回起こるまでの時間間隔を表す。機械が複数の箇所が壊れてから故障する場合や、指数分布に従って起こるイベントが複数回起こる場合に用いられる。無記憶性があるのは$k = 1$の時だけ。
2. 定義、性質
指数分布
$x > 0$のみで考える。単位時間あたりのイベント回数を$\lambda( > 0)$として、 $$ \text{pdf}(x; \lambda) = \lambda e^{-\lambda x} $$ $$ \text{cdf}(x; \lambda) = 1 - e^{-\lambda x} $$
あるいは平均間隔を$\beta (= 1/\lambda)$として、 $$ \text{pdf}(x; \beta) = \dfrac{1}{\beta} e^{-\frac{1}{\beta} x} $$ と表す人もいる。
期待値$= \dfrac{1}{\lambda} (= \beta)$、分散$ = \dfrac{1}{\lambda^{2}} (= \beta^{2})$
ガンマ分布*2
$x, k, \theta > 0$として、 $$ \text{pdf}(x; k, \theta) = \dfrac{\theta^{k} x^{k -1}e^{-\theta x}}{\Gamma(k)} $$ $$ \text{cdf}(x; k, \theta) = \dfrac{\gamma(k, \theta x)}{\Gamma(k)} $$
あるいは$\alpha = k$、$\beta = \dfrac{1}{\theta}$として、
$$ \text{pdf}(x; \alpha, \beta) = \dfrac{x^{\alpha-1}e^{-\frac{x}{\beta}}}{\beta^{\alpha}\Gamma(\alpha)} $$
ここで$\gamma(a, x)$は不完全ガンマ関数で$\gamma(a, x) = \displaystyle \int_{0}^{x} t^{a-1}e^{-t} dt$。$\theta$と$\beta$は定数倍しても$X \to cX$になるだけなので、scale parameter、一方$k$と$\alpha$はshape parameterと呼ばれる。
3. 導出
指数分布
- ポアッソン分布からの導出
ポアッソン分布は単位時間あたり$\lambda$回発生する事象が$k$回起こる確率を示し、 $$ \text{pdf}(X = k) = \dfrac{\lambda^{k}e^{-\lambda}}{\Gamma(k + 1)} $$
まだ事象が起こっていない時を考えるので$k \to 0$、さらに$x$までの時間間隔を考えたいので、$\lambda \to \lambda x$として、$e^{-\lambda x}$。時間$x$でイベントが起こる確率を知りたいので微分して符号を入れ替えると、$\lambda e^{-\lambda x}$が得られる。
- 幾何分布からの導出
幾何分布はベルヌーイ試行を繰り返して初めて成功するまでの回数Xの離散確率分布で、 $$ \text{pdf}(X = k) = p(1-p)^{k} $$
この分布を連続化するために離散的な試行を狭い範囲に無限に押しこめて、その代わりに一回あたりの失敗率も薄めてしまうことを考える。
ここで、$p \to \lambda \Delta_{x}$、$k \to \dfrac{x}{\Delta_{x}}$と置き換えことにしてそれぞれの変数の意味を考える。すると、$k$回で$kp = \lambda$回程度失敗するイベントが、$x$までの間の$\Delta_{x}$の区間$k$回では起こらず、$x + \Delta_{x}$で起こる確率を表すと解釈できる。知りたいのは確率密度関数なので$\Delta_{x}$で割って、
$$ \lim_{\Delta_{x}, p \to 0} \lambda(1 - p)^{\frac{\lambda x}{p}} = \lambda e^{-\lambda x} $$
参考:
幾何分布からの指数分布の導出: 指数分布の定義に間違いがあることに注意
指数分布 - NtRand
ガンマ分布
指数分布を畳み込み積分する。
参考: ときわ台学/統計学/ガンマ分布と指数分布
4. ガンマ分布が指数分布の共役事前分布であることを示す
ガンマ分布に従って発生させたパラメータは現時点までに複数回のイベントが起こるまでの時間を表すのであったが、ここではその意味を考えないことにする(としないと今のところ自分は理解できていない。)
事前分布は
尤度は $$ \text{Exp}(x_{1}, x_{2} \cdots x_{n}; \lambda) = \lambda^{n} e^{-\lambda \sum x_{i}} $$
なので、
$\alpha' = \alpha + n$、$\frac{1}{\beta'} = \frac{1}{\beta} + \sum x_{i}$と分かる。
参考: A Compendium of Conjugate Priorsの15ページ
5. プロット
指数分布
dexp(x, rate = 1, log = FALSE)となっていて、rateは$\lambda$に相当。
lambda <- c(0.5, 1, 1.5) for (i in 1:3) { plot(seq(0, 5, length=100), dexp(seq(0, 5, length=100), rate=lambda[i]), xlim=c(0,5), ylim=c(0, 1.5), type='l', col=c("red", "green", "blue")[i]) par(new=T) } # 最初の2引数は凡例の位置を指定 legend(3, 1.5, legend=c("lambda = 0.5", "lambda = 1", "lambda = 1.5"), lty=1, col=c("red", "green", "blue"))
ガンマ分布
dgamma(x, shape, rate = 1, scale = 1/rate, log = FALSE)となっていて、shapeが$k$、$\alpha$に、rateが$\beta$、scaleが$\theta$に相当。
k <- c(1, 2, 3, 5, 9) theta <- c(2, 2, 2, 1, 0.5) # 色はhttp://www.okada.jp.org/RWiki/?%BF%A7%B8%AB%CB%DCを参考にした for (i in 1:5) { plot(seq(0, 20, length=100), dgamma(seq(0, 20, length=100), shape=k[i], scale=theta[i]), xlim=c(0,20), ylim=c(0, 0.5), type='l', col=c("red", "green", "blue", "lightblue", "orange")[i]) par(new=T) } lgd <- rep("", 5) for (i in 1:5) { lgd[i] <- paste("k = ", as.character(k[i]),", θ = ", as.character(theta[i])) } legend(12, 0.5, legend=lgd, lty=1, col=c("red", "green", "blue", "lightblue", "orange"))
