背景
教科書が余りにも博物学的に分布を取り扱っている。この分布間には豊かな数学的構造や良い直観が含まれているのに......。
あと、偶に確率密度マンダラみたいなものを書いている人もいるけど、形式的な暗記に頼るのも良くない！
なので有名な分布について様々な観点で見つめ直してみる。そんな雑記。
1. まずは基本的な分布から
成功／失敗、表／裏、正常／不良のように2つの状態をとる試行を考える。一方のイベントが起こる確率を
$p$ とすれば、
f(k)=p^k(1-p)^{1-k},\qquad k=0,1
と表せる。これをベルヌーイ分布と呼ぶ。
「だから何やねん」と思うかもしれないが、$k=0, 1$のように状態を0/1に落とし込めるようになるというところが大事。
成功／失敗を ${0,1}$ に変換し、ベルヌーイ試行を $n$
回実施して得られる数列を $T_i$ とすると、
S_n=\sum_{i=1}^{n}T_i
が成功数になる。
このように総和を取れるようになるからエラいのだ。だから人類は確率変数の和を求めることができるかについて関心を払っている。
$n$ 回のうち $x$ 回成功する確率は、二項係数
_nC_x = \binom{n}{x}
を使って
P(X=x)=\binom{n}{x}p^x(1-p)^{n-x}.
これを二項分布と呼び、確率変数$X$が二項分布に従うことを
X\sim\operatorname{Bin}(n,p)
と書く。
1.1 r回成功するまでの試行回数
$T_r$ を $r$ 回成功するまでの試行回数とする。この時、
T_r\le n
とは「$n$ 回試行するまでに少なくとも $r$ 回は成功している」ことなので、
P(T_r\le n)=P(S_n\ge r).
$S_n\sim\operatorname{Bin}(n,p)$ より、
P(T_r\le n)=\sum_{k=r}^{n}\binom{n}{k}p^k(1-p)^{n-k}.
知りたいのは $P(T_r=n)$ なので、
P(T_r=n)=P(T_r\le n)-P(T_r\le n-1).
パスカルの関係
\binom{n}{k}=\binom{n-1}{k-1}+\binom{n-1}{k}
を使って整理すると、
P(T_r\le n)
=\binom{n-1}{r-1}p^r(1-p)^{n-r}
+P(T_r\le n-1).
したがって、
P(T_r=n)=\binom{n-1}{r-1}p^r(1-p)^{n-r}.
この確率関数を負の二項分布と呼び、
T_r\sim\operatorname{NB}(r,p)
と書く。特に $r=1$ のとき、
P(T_1=n)=p(1-p)^{n-1}
となり、これを幾何分布と呼び、
T_1\sim\operatorname{Geo}(p)
と書く。
「$n$ 回の内 $x$ 回が成功」と「$r$ 回成功するために必要な試行回数」は、互いにひっくり返った関係にある。このような関係をここでは双対性と呼ぶ。
2. もっと細かくしたい！：二項分布からポアソン分布へ
いままでは、試行回数や成功回数に着目してきた訳だが、例えば、
「単位時間当りにイベントが起こる確率」
を考えることは可能だろうか？
ナイーブな考え方として、「100回試行する」ということを「100秒間、毎秒観測する」と捉えれば穏当そうだ（つまり毎秒試行している）。ただここで問題がある。
いままでのモデルでは、試行1回につき $0/1$ の一方のみが1つ返ってくることが要求されている。つまり、1秒間に複数のイベントが起こってはいけないのだ！
「じゃあもっと単位時間を細かく（例えば $1/1000$ 秒とか）すればいいじゃん？」
となる訳だが、そうすると単位時間当りのイベント発生確率 $p$ も十分小さくならないとベルヌーイ試行を近似できない。
例えば、毎秒 $0.1$ の確率で発生するイベントを $1/1000$ 秒単位で見るなら、その細かさに合わせて1区間あたりの確率も小さくする必要がある。
この場合、$0.1 \times \dfrac{1}{1000} = \dfrac{1}{10000}$となるであろう。ここで試行回数と生起確率の積は一定で維持されている。
これを定式化しよう。
単位時間を $\Delta t$、発生率を $\lambda$ とし、
p=\lambda\Delta t
とする。単位時間当りの試行回数を
n=\frac{1}{\Delta t}
とすれば、
np=\frac{1}{\Delta t}\lambda\Delta t=\lambda
で、$\Delta t$ の細かさに依らず一定である。
長さ $T$ の時間を観測すると、時間区間の個数は $T/\Delta t$ なので、
X_{\Delta t}\sim\operatorname{Bin}\left(\frac{T}{\Delta t},\lambda\Delta t\right).
したがって、
P(X_{\Delta t}=x)
=\binom{T/\Delta t}{x}(\lambda\Delta t)^x
(1-\lambda\Delta t)^{T/\Delta t-x}.
$\Delta t\to0$ を考える。組合せ部分は
\binom{T/\Delta t}{x}(\lambda\Delta t)^x
\to \frac{(\lambda T)^x}{x!},
また、
(1-\lambda\Delta t)^{T/\Delta t}\to e^{-\lambda T},
\qquad
(1-\lambda\Delta t)^{-x}\to1.
よって、
P(X=x)=\frac{(\lambda T)^x}{x!}e^{-\lambda T}.
これをポアソン分布と呼ぶ。
2.1 負の二項分布を連続化すると？
負の二項分布
P(N_r=n)=\binom{n-1}{r-1}p^r(1-p)^{n-r}
に $p=\lambda\Delta t$ を代入すると、
P(N_r=n)=\binom{n-1}{r-1}(\lambda\Delta t)^r
(1-\lambda\Delta t)^{n-r}.
時刻 $t$ に対応する試行番号を
n=\frac{t}{\Delta t}
とする。$\Delta t\to0$ では、これは $r$ 回目の成功までにかかる総待ち時間に相当する。時間は連続変数なので、
f(t)=\lim_{\Delta t\to0}\frac{1}{\Delta t}
P\left(N_r=\frac{t}{\Delta t}\right)
と密度を定義する。
まず、
\binom{n-1}{r-1}
=\frac{(n-1)(n-2)\cdots(n-r+1)}{(r-1)!}.
$n\Delta t=t$ より、
\binom{n-1}{r-1}(\Delta t)^{r-1}
=
\frac{(t-\Delta t)(t-2\Delta t)\cdots(t-(r-1)\Delta t)}
{(r-1)!}
\to
\frac{t^{r-1}}{(r-1)!}.
また、
(1-\lambda\Delta t)^{n-r}
=(1-\lambda\Delta t)^{t/\Delta t}(1-\lambda\Delta t)^{-r}
\to e^{-\lambda t}.
したがって、
f(t)
=\frac{t^{r-1}}{(r-1)!}\lambda^r e^{-\lambda t}.
これをアーラン分布と呼び、$\operatorname{Erlang}(r,\lambda)$ と書く。
ここで、アーラン分布では $r$ は「$r$ 回目のイベントまで待つ」という意味なので整数である。
一方、密度関数
f(t)
=
\frac{t^{r-1}}{(r-1)!}\lambda^r e^{-\lambda t}
だけを見ると、$(r-1)!$ をうまく一般化できれば、整数だった $r$ を連続的なパラメータへ拡張できそうである。
そこで、階乗を一般化するガンマ関数
\Gamma(\alpha)
=
\int_0^\infty x^{\alpha-1}e^{-x}\,dx
を導入する。
正の整数 $r$ に対して、
\Gamma(r)=(r-1)!
が成り立つので、アーラン分布は
f(t)
=
\frac{t^{r-1}}{\Gamma(r)}
\lambda^r e^{-\lambda t}
とも書ける。
ここで整数 $r$ を正の実数 $\alpha$ に拡張すると、
f(t)
=
\frac{t^{\alpha-1}}{\Gamma(\alpha)}
\lambda^\alpha e^{-\lambda t}
となる。
これをガンマ分布と呼び、$\operatorname{Ga}(\alpha,\lambda)$ と書く。
Column：ガンマ関数とは何者か？
ガンマ関数は
\Gamma(\alpha)
=
\int_0^\infty x^{\alpha-1}e^{-x}\,dx
と定義される。
\Gamma(\alpha+1)
=
\int_0^\infty x^\alpha e^{-x}\,dx
を部分積分すると、
\Gamma(\alpha+1)
=
\left[-x^\alpha e^{-x}\right]_0^\infty
+
\alpha\int_0^\infty x^{\alpha-1}e^{-x}\,dx.
境界項は $0$ なので、
\Gamma(\alpha+1)=\alpha\Gamma(\alpha).
また、
\Gamma(1)=\int_0^\infty e^{-x}\,dx=1.
従って、
\Gamma(2)=1\Gamma(1)=1!,
\Gamma(3)=2\Gamma(2)=2!,
と続き、一般に
\Gamma(n)=(n-1)!
となる。アーラン分布の $(r-1)!$ を $\Gamma(r)$
に置き換えることで、整数 $r$ から実数 $\alpha$
への拡張が自然に見える。
ここまでをまとめると、
\text{ベルヌーイ分布}
\xrightarrow{\text{和を取る}}
\text{二項分布}
\xrightarrow{\text{連続化}}
\text{ポアソン分布},
および
\text{二項分布}
\xrightarrow{\text{双対}}
\text{負の二項分布}
\xrightarrow{\text{連続化}}
\text{アーラン分布}
\longrightarrow
\text{ガンマ分布}
という構造が見える。
2.2 ポアソン分布とガンマ分布の関係
このように見てみるとかなりきれいな構造である。そこで、ポアソン分布とガンマ分布の関係性についても眺めてみよう。
ポアソン分布は
P_0(x)
=\frac{(\lambda t)^x}{x!}e^{-\lambda t}
=\frac{(\lambda t)^x}{\Gamma(x+1)}e^{-\lambda t},
ガンマ分布は
\operatorname{Ga}(\alpha,\lambda)
=\frac{t^{\alpha-1}}{\Gamma(\alpha)}
\lambda^\alpha e^{-\lambda t}
で、よく似た式である。
$N(t)$ を時刻 $t$ までのイベント回数、$T_r$ を $r$
回目のイベントが発生する時刻と置くと、
T_r\le t
\iff
N(t)\ge r.
つまり、
P(T_r\le t)=P(N(t)\ge r).
右辺はポアソン分布なので、
F_{T_r}(t)
=1-\sum_{k=0}^{r-1}\frac{(\lambda t)^k}{k!}e^{-\lambda t}.
すなわち、
F_{T_r}(t)
=
1-e^{-\lambda t}
\left(
\frac{(\lambda t)^0}{0!}
+\frac{(\lambda t)^1}{1!}
+\cdots
+\frac{(\lambda t)^{r-1}}{(r-1)!}
\right).
$t$ で微分すると、中間項が打ち消し合って、
f_{T_r}(t)
=
\frac{\lambda^r t^{r-1}}{(r-1)!}e^{-\lambda t}
=
\frac{\lambda^r t^{r-1}}{\Gamma(r)}e^{-\lambda t}.
したがって、
T_r\sim\operatorname{Ga}(r,\lambda)
が得られる。
これは負の二項分布を求めた時の流れとまったく同じである（むしろ微分を使える分、計算はラク）。
つまり、ポアソン分布とガンマ分布もまた双対の関係だったのである。
3. 順番も大事！：一様分布からベータ分布へ
ここまでは「イベントが何回発生したか？」に注目していた。しかし実際には、あるしきい値を基に成功／失敗を判断するケースも多い。
そこで、発生したものを大きい順あるいは小さい順に並べられることを考える。
簡単のために $[0,1]$ の区間で生成される一様乱数を考える。この乱数機で$n$ 個の乱数
U_1,U_2,\ldots,U_n
が得られたとする。
$x\in[0,1]$ をしきい値としたとき、$x$ 以下である乱数の個数 $N(x)$ は
N(x)\sim\operatorname{Bin}(n,x)
に従う。
この乱数列を小さい順に並べかえた時に得られる数列を
U_{(1)},U_{(2)},\ldots,U_{(n)}
とする。
U_{(1)}\le U_{(2)}\le\cdots\le U_{(n)}.
ここで $U_{(r)}$ を $r$ 番目の順序統計量と呼ぶ。
$U_{(r)}\le x$ を考える。これは「$x$ 以下の乱数が少なくとも $r$
個ある」ことと同値なので、
U_{(r)}\le x
\iff
N(x)\ge r.
したがって、
P(U_{(r)}\le x)=P(N(x)\ge r).
$N(x)\sim\operatorname{Bin}(n,x)$ より、
F_{U_{(r)}}(x)
=
\sum_{k=r}^{n}\binom{n}{k}x^k(1-x)^{n-k}.
これを $x$ で微分すれば $P(U_{(r)}=x)$ に対応する密度が得られる。
数学者が近くに居ないことを確認し、いきなりシグマの中身を微分してしまおう（問題が起こらなければ良いのである）。
\begin{aligned}
f_{U_{(r)}}(x)
&=
\sum_{k=r}^{n}
\left[
k\binom{n}{k}x^{k-1}(1-x)^{n-k}
-(n-k)\binom{n}{k}x^k(1-x)^{n-k-1}
\right].
\end{aligned}
ここで、
k\binom{n}{k}=n\binom{n-1}{k-1},
\qquad
(n-k)\binom{n}{k}=n\binom{n-1}{k}
を使うと、
f_{U_{(r)}}(x)
=
n\sum_{k=r}^{n}
\left[
\binom{n-1}{k-1}x^{k-1}(1-x)^{n-k}
-
\binom{n-1}{k}x^k(1-x)^{n-k-1}
\right].
隣接項が順々に打ち消し合い、端の項だけが残るので、
f_{U_{(r)}}(x)
=
n\binom{n-1}{r-1}x^{r-1}(1-x)^{n-r}.
すなわち、
f_{U_{(r)}}(x)
=
\frac{n!}{(r-1)!(n-r)!}
x^{r-1}(1-x)^{n-r}.
ガンマ関数を使えば、
f_{U_{(r)}}(x)
=
\frac{\Gamma(n+1)}
{\Gamma(r)\Gamma(n-r+1)}
x^{r-1}(1-x)^{n-r}.
ベータ関数
B(a,b)=\frac{\Gamma(a)\Gamma(b)}{\Gamma(a+b)}
を使うと、通常のパラメータ表示では
f_{U_{(r)}}(x)
=
\frac{x^{r-1}(1-x)^{n-r}}
{B(r,n-r+1)},
\qquad 0<x<1.
したがって、
U_{(r)}\sim\operatorname{Beta}(r,n-r+1).
つまり、一様分布から得た標本を並べ、その $r$
番目を見るとベータ分布が現れる。
Column：ベータ関数とは何だったのか
ベータ分布を導出したところで、突然
B(a,b)=\frac{\Gamma(a)\Gamma(b)}{\Gamma(a+b)}
というベータ関数が出てきたので、少し補足しておく。
ベータ関数は本来、
B(a,b)
=
\int_0^1 x^{a-1}(1-x)^{b-1}\,dx,
\qquad a,b>0
と定義される関数である。
つまり、
x^{a-1}(1-x)^{b-1}
を $0$ から $1$ まで積分した値そのものに名前を付けたものだ。
なぜこんなものにわざわざ名前を付けるのか？
ベータ分布の形
x^{a-1}(1-x)^{b-1}
をそのまま確率密度関数にしようとしても、一般には積分値が $1$
にならない。
そこで、その全積分
\int_0^1 x^{a-1}(1-x)^{b-1}\,dx
=
B(a,b)
で割ってやれば、
\int_0^1
\frac{x^{a-1}(1-x)^{b-1}}{B(a,b)}
\,dx
=1
となる。
つまりベータ分布における $B(a,b)$ は、ひとことで言えば
$x^{a-1}(1-x)^{b-1}$ を確率密度関数にするための正規化定数
である。
Column：ガンマ分布とベータ分布もつながっている
ここまで、ガンマ分布とベータ分布は別々の経路から登場した。
しかし実は、この2つの分布にはかなり直接的な関係がある。
独立な確率変数
X\sim\operatorname{Ga}(\alpha,\lambda),
\qquad
Y\sim\operatorname{Ga}(\beta,\lambda)
を考える。ただし、2つのガンマ分布の rate $\lambda$ は同じものとする。
ここから、
S=X+Y,
\qquad
R=\frac{X}{X+Y}
という2つの新しい確率変数を作ってみる。
$S$ は $X,Y$ を合わせた全体量、$R$ はその全体のうち $X$
が占める割合である。
逆に書けば、
X=RS,
\qquad
Y=(1-R)S.
つまり、$(X,Y)$ という「2つの量」を、$(S,R)$
という「全体量と割合」に座標変換している。
この変数変換のヤコビアンを計算すると、
\frac{\partial(X,Y)}{\partial(S,R)}
=
\begin{vmatrix}
\dfrac{\partial X}{\partial S} & \dfrac{\partial X}{\partial R}\\[6pt]
\dfrac{\partial Y}{\partial S} & \dfrac{\partial Y}{\partial R}
\end{vmatrix}
=
\begin{vmatrix}
R & S\\
1-R & -S
\end{vmatrix}
=-S.
したがって、その絶対値は
\left|
\frac{\partial(X,Y)}{\partial(S,R)}
\right|
=S.
$X,Y$ は独立なので、同時密度は
f_{X,Y}(x,y)
=
\frac{\lambda^\alpha}{\Gamma(\alpha)}
x^{\alpha-1}e^{-\lambda x}
\frac{\lambda^\beta}{\Gamma(\beta)}
y^{\beta-1}e^{-\lambda y}.
ここに
x=rs,
\qquad
y=(1-r)s
とヤコビアン $s$ を代入すると、
\begin{aligned}
f_{S,R}(s,r)
&=
f_{X,Y}(rs,(1-r)s)
\left|
\frac{\partial(X,Y)}{\partial(S,R)}
\right|
\\
&=
\frac{\lambda^{\alpha+\beta}}
{\Gamma(\alpha)\Gamma(\beta)}
(rs)^{\alpha-1}
\{(1-r)s\}^{\beta-1}
e^{-\lambda s}s
\\
&=
\frac{\lambda^{\alpha+\beta}}
{\Gamma(\alpha)\Gamma(\beta)}
s^{\alpha+\beta-1}e^{-\lambda s}
r^{\alpha-1}(1-r)^{\beta-1}.
\end{aligned}
ここで、$\Gamma(\alpha+\beta)$ を掛けて割ると、
\begin{aligned}
f_{S,R}(s,r)
&=
\underbrace{
\frac{\lambda^{\alpha+\beta}}
{\Gamma(\alpha+\beta)}
s^{\alpha+\beta-1}e^{-\lambda s}
}_{\operatorname{Ga}(\alpha+\beta,\lambda)}
\\
&\qquad\times
\underbrace{
\frac{\Gamma(\alpha+\beta)}
{\Gamma(\alpha)\Gamma(\beta)}
r^{\alpha-1}(1-r)^{\beta-1}
}_{\operatorname{Beta}(\alpha,\beta)}.
\end{aligned}
きれいに「$S$ だけの関数」と「$R$ だけの関数」に分離した。
したがって、
\boxed{
S=X+Y\sim\operatorname{Ga}(\alpha+\beta,\lambda)
}
かつ、
\boxed{
R=\frac{X}{X+Y}\sim\operatorname{Beta}(\alpha,\beta)
}
となる。
さらに同時密度が
f_{S,R}(s,r)=f_S(s)f_R(r)
と積に分離しているので、
\boxed{S\perp R}
でもある。
つまり、「全体がどれだけ大きいか」と「そのうち $X$
が何割を占めるか」は独立なのである。
また、ベータ分布側の正規化定数を見れば、
\frac{1}{B(\alpha,\beta)}
=
\frac{\Gamma(\alpha+\beta)}
{\Gamma(\alpha)\Gamma(\beta)}
でなければならない。従って、
\boxed{
B(\alpha,\beta)
=
\frac{\Gamma(\alpha)\Gamma(\beta)}
{\Gamma(\alpha+\beta)}
}
という、先ほど使ったベータ関数とガンマ関数の関係もここで同時に得られる。
さらに今回、順序統計量から
U_{(r)}
\sim
\operatorname{Beta}(r,n-r+1)
を得ていた。
一方、
X\sim\operatorname{Ga}(r,\lambda),
\qquad
Y\sim\operatorname{Ga}(n-r+1,\lambda)
とすれば、
\frac{X}{X+Y}
\sim
\operatorname{Beta}(r,n-r+1).
したがって分布の意味で、
U_{(r)}
\overset{d}{=}
\frac{X}{X+Y}
と書ける。
一様分布の順序統計量から現れたベータ分布と、待ち時間から現れたガンマ分布が、ここで合流する。