グレた確率統計 ～ポアソン分布～

これまでにしてきたことを少し振り返ってみましょう。

まず、ベルヌーイ試行から始まりました。ONかOFFか、YESかNOか、1か0か、というやつですね。それを複数回行ったときのONの数の確率分布を求めました（二項分布）。そして、初めてONが出たときの試行回数の確率分布も求めました（幾何分布）。畳み込みを踏まえた上で、幾何分布の一般形、すなわち r回ONが出たときの試行回数の確率分布も求めました（負の二項分布）。そのあと、連続化に進み、幾何分布を連続化しました（指数分布）。また、負の二項分布の連続化、もしくは指数分布の畳み込みを行いました（アーラン分布）。

二項分布、幾何分布、負の二項分布、指数分布、アーラン分布が、単純なベルヌーイ試行から導出されてきたことを見てきたわけです。

さて、今回は、次のような分布を求めてみます。すなわち、「単位時間で事象がk回起こる確率」です。

最も単純明快な導出方法は、おそらく、アーラン分布（ガンマ分布）を用いたものでしょう。

単位時間で事象がk回起こる確率というのは、t' < 1 なる t' でk回起こり、残りの時間は起こらなかった確率と言い換えられます。t = t' で k回起こる確率というのはガンマ分布、

$$\Gamma(t; r=k) = \frac{t^{k-1} \phi^{k} \exp(-\phi t)}{(k-1)!}$$

で表わせ、残りの時間、すなわち t = 1-t' で1回も起こらない確率$g_0(t)$は指数分布を用いて、

$$g_0(t) = 1-\int_0^t \phi \exp(-\phi t') dt' = \exp(-\phi t)$$

とできますから、この2つの畳込み、

$$Pr(X=k) = \int_0^1 \Gamma(t';r=k)\cdot g_0 (1-t')dt' \\ = \int_0^1 \frac{t'^{k-1} \phi^{k} \exp(-\phi t')}{(k-1)!} \cdot \exp(-\phi (1-t')) dt' \\ = \frac{\phi^k \exp(-\phi)}{(k-1)!} \int_0^1 t'^{k-1}dt' = \frac{\phi^k}{k!}  \exp(-\phi) $$

この分布をポアソン分布といい、意味としてはさきに書いた通り、「単位時間で事象がk回起こる確率」となります。

というわけで、主要な確率分布が以上のようにすべて繋がっているように考えられるということが分かります。

（多分おしまい）

グレた確率統計 ～ガンマ・アーラン分布～

前回は、幾何分布から指数分布へ、つまり離散型から連続型への変換をやりました。
幾何分布は畳み込むことで負の二項分布になりましたね。
もちろんのことながら、負の二項分布を連続化した分布も考えることができます。
さらに言えば、指数分布を畳み込むことで同じ分布が得られるはずです。
それこそがガンマ分布（厳密に言えばアーラン分布）なわけです。

そこで、ここではガンマ分布を2つのやり方
① 負の二項分布の連続化
② 指数分布の畳み込み
で導出してみます。

① 負の二項分布の連続化

負の二項分布は、
$$Pr(X=k) = ~{}_{k-1} \mathrm{C}_{r-1} ~ p^r q^{k-r}$$
と表わせ、「事象が試行k回目にちょうどr回生じる確率」の意味付けができます。
これを連続化すれば、おそらく
「事象が時間tのときにちょうどr回生じる確率分布」が得られるはずです。

まず、組み合わせの変形はそのままではやりづらいので、負の二項分布の自然対数を変形していくことにします。
また、非常に大きい数 x について、スターリンの近似式、
$$ \ln x! \approx x \ln x - x$$
が成り立つこと、そして、小さい数 x について、次の近似式が成り立ちます：
$$ \ln (1-x) \approx -x $$
$$\ln Pr(X=k) = \ln \left( ~{}_{k-1} \mathrm{C}_{r-1} ~ p^r q^{k-r} \right) = \ln \left( \frac{(k-1)!}{(k-r)! (r-1)!} ~ p^r q^{k-r} \right ) \\ = \ln{(k-1)!} ~- \ln(k-r)! ~-\ln(r-1)! ~+ r\ln p + (k-r)\ln q $$

スターリンの式で、$\ln(k-1)!$と$\ln(k-r)!$を近似し、$p=\phi dt$を代入すると

$$= (k-1)\ln(k-1)-(k-1) - (k-r)\ln(k-r) +(k-r) -\ln(r-1)!+r \ln(\phi dt)+(k-r)\ln(1-p)$$

ここで、次の近似、
$$\ln(k-a)=\ln k + \ln(1-a/k) \approx \ln k - \frac{a}{k} \\\ln(1-p) \approx -p = -\phi dt$$
を使うと、
$$= (k-1)\left(\ln k-\frac{1}{k}\right) - (k-r)\left(\ln k -\frac{r}{k}\right) +1-r -\ln(r-1)!+r \ln(\phi dt)-(k-r)\phi dt\\= (r-1)\ln k -\frac{1}{k} ((k-1)-r(k-r)) + 1-r -\ln(r-1)!+ r\ln\phi +r\ln dt -\phi t+r\phi dt$$
少し式が煩雑になってきたので、
$L = -\ln(r-1)!+r\ln\phi-\phi t$と置いて、式を整理していくと、
$$=(r-1)\ln k -\frac{1}{k} ((k-1)-r(k-r)) + 1-r +r\ln dt+r\phi dt + L\\=(r-1)(\ln t - \ln dt)-\frac{1+r^2}{k}+r\ln dt + r\phi dt +L\\=(r-1)\ln t + \ln dt - \left(\frac{(1+r)^2}{t} +r\phi \right )dt + L$$
ここで、Lというのは、
$$L = \ln \left( \frac{\phi^r \exp(-\phi t)}{(r-1)!}  \right )$$
と変形できるので、これを使って変形すると、
$$= \ln \left( \frac{t^{r-1} \phi^r \exp(-\phi t)dt}{(r-1)!} \right) - \left(\frac{(1+r)^2}{t} +r\phi \right )dt\\ \approx \ln \left( \frac{t^{r-1} \phi^r \exp(-\phi t)dt}{(r-1)!} \right) ~~(\because dt \ll 1)$$

よって、
$$Pr(T=t) = \frac{t^{r-1} \phi^r \exp(-\phi t)}{(r-1)!} dt$$
となり、$Pr(T=t) = \Gamma (t) dt $とすれば、確率密度 Γ(t) は、
$$ \Gamma (t) = \frac{t^{r-1} \phi^r \exp(-\phi t)}{(r-1)!} $$
となります。

② 指数分布の畳み込み

t秒後にちょうどk回起こる確率分布を$g_k(t)$とすると、k=2のときは2つの指数分布をたたみ込めばいいので、

$$g_2(t)=\int_0^t \left(\phi \exp(-\phi (t-t')) \right )\left(\phi \exp(-\phi t') \right )dt'\\ = \phi^2 \exp(-\phi t) \int_0^t ~dt' = t \phi^2 \exp(-\phi t)$$

となります。一応$g_3(t)$も求めてみると、これは$g_2(t)$と指数分布の畳み込みで求められ、

$$ g_3(t) = \int_0^t g_2(t') \phi \exp(-\phi(t-t')) dt' \\ = \int_0^t t' \phi^3 \exp(-\phi t) dt' = \frac{t^2 \phi^3 \exp(-\phi t)}{2!} $$

となります。

ここから帰納的に、

$$ g_k(t) = \frac{t^{k-1} \phi^k \exp(-\phi t)}{(k-1)!} $$

と仮定してみると、

$$g_{k+1}(t) = \int_0^t g_k(t') \phi \exp(-\phi(t-t'))dt' \\ = \int_0^t \frac{t'^{k-1} \phi^k \exp(-\phi t')}{(k-1)!} \cdot \phi \exp(-\phi(t-t')) dt' \\ = \frac{\phi^{k+1} \exp(-\phi t)}{(k-1)!} \int_0^t t'^{k-1} dt' = \frac{\phi^{k+1} \exp(-\phi t)}{(k-1)!} \cdot \frac{t^k}{k}\\ = \frac{t^k \phi^{k+1} \exp(-\phi t)}{k!} $$

ということで、やはりどちらの側からもガンマ分布が求められました。
ちなみに今回出した分布のパラメータ r は正の整数値をとりますが、実際のガンマ分布のパラメータ r は正の実数値を取ります。

そこで、正の整数値をとる場合の確率分布はアーラン分布と呼ばれています。ですから、これら一連の導出は厳密にはアーラン分布の導出ということになります。

さて、ガンマ分布の平均値と分散を求めてみましょう。愚直に計算してもいいのですが、どうせなら、負の二項分布を連続化、もしくは指数分布を畳み込んだものがガンマ分布ということを利用してみましょう。

負の二項分布の平均値と分散はそれぞれ、$\mu = r/p,~~ \sigma^2 = rq/p^2$となります。
ここに、変数変換 $T = X dt$を考慮しながら、 $p = \phi dt$ を代入すれば、
$$ E(T) = E(X dt) = E(X) dt = \frac{r dt}{p} = \frac{r}{\phi} \\ Var(T) = Var(X dt) = Var(X) dt^2 = \frac{r q dt^2}{p^2} = \frac{r}{\phi^2} $$
と求められます。

また、ガンマ分布に従う確率変数 Z と r個の指数分布に従う変数 $X_i$ は、
$ Z = X_1 + X_2 + ... + X_r$
であり、かつ、$X_i$ はそれぞれ独立なので、
$E(Z) = E(X_1 + X_2 + ... + X_r) = E(X_1) + E(X_2)+ ... + E(X_r) = r/\phi$
分散も同様に、
$Var(Z) = Var(X_1 + X_2 + ... + X_r) = Var(X_1) + Var(X_2)+ ... + Var(X_r) = r/\phi^2$
と求められます。

定義通りに求めるやり方は、練習問題としておきます。

グレた確率統計 ～指数分布の補足～

さて、前回の幾何分布と指数分布について少し補足したいと思います。


指数分布の導出過程で、
$$p = \frac{dt}{\mu}$$
と置きましたね。すなわち、ここでpは「微小時間間隔 dt で事象が起こる確率」となっています。
これは、T = X dt と置いたことからも分かります。
この式を少し眺めていると、

[確率] = [確率密度] × [微小時間]

の形に見えてきませんか？つまり、pと対応する確率密度 $\phi$ を、

$$ p = \phi dt $$
$$ \phi = p / dt = 1/\mu $$

と定義してやることができます。これを指数分布に代入してみると、
$$f(t) = \frac{1}{\mu}\cdot \mathrm{exp}(-t/\mu) = \phi \cdot \exp(-\phi t)$$
となります。指数分布の離散版ともいえる幾何分布が、
$$Pr(X=k)= pq^{k-1} ~~~ (p+q=1)$$
だったことを考えると、qの累乗と、指数部に対応付けがしたくなります。

具体例を考えてみると、たとえば幾何分布では、
$$Pr(X \geq  k) = 1 - \sum_{i=1}^{k}p(1-p)^{i-1} \\ = 1 - p\cdot\frac{1-q^{k}}{1-q} = 1- (1-q^k) = q^k$$
であり、これに対応する指数分布の場合を考えると、
$$Pr(T \geq  t) = 1 - \int_{0}^{t} \phi \exp(-\phi t') dt'\\= 1 - [-\exp(-\phi t')]_0^t = 1- (1- \exp(-\phi t)) = \exp(-\phi t)$$
となり、
$$q^k = \exp(-\phi t)$$
と関連づけれそうです。
導出過程のことを考えると、tにとって k と k-1 はほとんど変わらないので、
$$ f(t) = \phi \cdot \exp(-\phi t) $$
というのは、幾何分布において、
$$p \rightarrow \phi$$
$$q^{k-1} \rightarrow \exp(-\phi t)$$
と置き換えたものだとみれます。

もうひとつ、q と exp(-φt) の対応付けが妥当な例を見てみましょう：
$$ q^{k-1} = (1-p)^{k-1} = (1-\phi dt)^k = (1-\phi dt)^{t/dt} = \exp(-\phi t)$$
ここで、$1/({1-\phi dt}) \approx 1$, $k = t/dt$ を用いました。
最後の等式は、指数分布の導出でも用いた、ネイピア数の定義です。

ということで、幾何分布と指数分布の形式的な対応付けが与えられたことになります。
その結果、期待値と分散も幾何分布から指数分布への簡単な対応付けで得られるはずです。

グレた確率統計 ～指数分布～

今までは離散分布だけを扱ってきました。
ここでひとつ、幾何分布を連続化させてみましょう。
というのも、幾何分布の意味するところは「何かが起こるまでにかかる回数」なわけですが、
これを「何かが起こるまでにかかる時間」とするのはかなり有意義そうだからです。

まず、幾何分布は、
$$ Pr(X=k) = p q^{k-1}$$
で、期待値と分散は、
$$ E[X] = \frac{1}{p}$$
$$ Var[X] = \frac{q}{p^2}$$
となりましたね。

ここで、時間確率変数 T を微小量 $\Delta t$ を用いて、$T = \Delta t \cdot X$ と定義しましょう。
そうすると、幾何分布はtと $\Delta t$ で書き換えられて、

$$ Pr(T=t)= Pr(\Delta t \cdot X = t) = Pr(X = t/\Delta t) = pq^{(t/\Delta t) -1} $$

となります。ここで、Tの期待値を求めてみると、
$$\mu = E[T]=E[X\Delta t]=\Delta t \cdot E[X]=\frac{\Delta t}{p} $$
となります。

次に、p を $\mu$ で表すと、
$$p = \frac{\Delta t}{\mu}$$
となります。これを先ほどの分布の式に代入してやると、
$$ Pr(T=t)= \frac{\Delta t}{\mu}~\left(1-\frac{\Delta t}{\mu}\right)^{(t/\Delta t) -1} $$

さて、連続分布なので、確率密度関数を
$$f(t) = \frac{Pr(T=t)}{\Delta t}$$
で定義すると、

$$ f(t) = \frac{1}{\mu}~\left(1-\frac{\Delta t}{\mu}\right)^{(t/\Delta t) -1}\\
= \frac{1}{\mu (1-(\Delta t/\mu))}\left( \left(1-\frac{\Delta t}{\mu} \right )^{-\mu/\Delta t}\right)^{-t/\mu} \\
\xrightarrow[]{\Delta t \rightarrow 0} \frac{1}{\mu}\cdot \mathrm{exp}(-t/\mu)$$

この最後の式、
$$f(t) = \frac{1}{\mu}\cdot \mathrm{exp}(-t/\mu)$$
は指数分布といい、もちろんのことながら、「何かが起こるまでにかかる時間」の確率分布です。

期待値と分散は（導出は練習問題としましょう！）
$$ E[X] = \mu $$
$$ Var[X] = \mu ^2 $$
さらに、指数分布も幾何分布と同様、無記憶性を持っています（これも練習問題としましょう）。

グレた確率統計 ～忌まわしき組み合わせ（補遺）～

次の記事から、いよいよ離散分布から連続分布へ移っていきます。
そのターニングポイントとして、ここで天下り的に用いた３つの組み合わせの公式（とその補題）を、
ここで示しておこうかと思います。
というのは、
1. 負の二項定理（？）
$$\sum_{r=0}^\infty {}_{n+r-1}\mathrm{C}_{r}~x^r = (1-x)^{-n}$$

2.
$${ }_{n+m}\mathrm{C}_k =\sum_{l = k-\mathrm{min}(k,m)}^{\mathrm{min}(k,n)} { }_n\mathrm{C}_l \cdot { }_m\mathrm{C}_{k-l}$$

3.
$${}_{m}\mathrm{C}_{m}+  {}_{m+1}\mathrm{C}_{m}+{}_{m+2}\mathrm{C}_{m}+...+{}_{n}\mathrm{C}_{m}= { }_{n+1}\mathrm{C}_{m+1}$$

そして、パスカルの三角形を数式で表した、
4.
$${}_{n}\mathrm{C}_{m} = { }_{n-1}\mathrm{C}_{m} + { }_{n-1}\mathrm{C}_{m-1}$$

を補題とし、計４題を示していきたいと思います。

# 1. の証明
これは、$\frac{1}{(1-z)^n}$を級数展開してみればよい。
まずは、導関数を求める。
$$\frac{\mathrm{d} }{\mathrm{d} x} (1-z)^{-n} = n(1-z)^{-n-1}\\
\frac{\mathrm{d^2} }{\mathrm{d} x^2} (1-z)^{-n} =\frac{\mathrm{d} }{\mathrm{d} x} n(1-z)^{-n-1} = n(n+1)(1-z)^{-n-2}\\$$
これより、帰納的に
$$\frac{\mathrm{d}^k }{\mathrm{d} x^k} (1-z)^{-n} = n(n+1)...(n+k-1)(1-z)^{-n-k}=\frac{(n+k-1)!}{(n-1)!}(1-z)^{-n-k}$$
と分かる。
これを用いてz = 0 のまわりで級数展開すると、
$$(1-z)^{-n} = \sum_{k=0}^\infty \frac{(n+k-1)!}{(n-1)!~k!} z^k = \sum_{k=0}^\infty { }_{n+k-1}\mathrm{C}_{k}\cdot z^k$$
ここで、この級数の収束半径を求める。第k項を$a_k$と書くとすると、
ダランベールの収束判定法を用いて、
$$\left| \frac{a_{k+1}}{a_k} \right| = \left| \frac{n+k}{k+1}~ z \right| \xrightarrow[]{k \rightarrow \infty} |z| $$
よって、収束半径は $|z| < 1$ である。zが確率であるならば、$0 \leq z \leq 1$であるから、自明な場合の多い z = 1 を除けば、この等式は常に成り立つ。

# 2. の証明
$$(1+x)^{n+m} = \sum_{k=0}^{n+m} { }_{n+m}\mathrm{C}_{k} \cdot x^k\\
(1+x)^n (1+x)^m = \left( \sum_{i=0}^{n} { }_{n}\mathrm{C}_{i} \cdot x^i \right ) \left( \sum_{j=0}^{m} { }_{m}\mathrm{C}_{j} \cdot x^j \right )$$

$(1+x)^{n+m} = (1+x)^n (1+x)^m$ であるから、$x^k$の係数は、
(左辺) $$= { }_{n+m}\mathrm{C}_{k}$$
(右辺) $$=\sum_{\substack{i+j = k \\ 0 \leq i \leq n \\ 0 \leq j \leq m}} { }_{n}\mathrm{C}_{i} \cdot { }_{m}\mathrm{C}_{j}$$
$0 \leq k \leq m+n$より、これはすなわち
$$= \sum_{i = \mathrm{max}(k-m,0)}^{\mathrm{min}(n,k)} { }_{n}\mathrm{C}_{i} \cdot { }_{m}\mathrm{C}_{k-i}$$
であるので、題意は満たされた。ちなみに、
$$k-\mathrm{min}(k,m) = k + \mathrm{max} (-k, -m) = \mathrm{max}(0, k-m) $$
となる。

# 4. の証明
まずは補題として，パスカルの三角形を証明します。

$${ }_{n-1}\mathrm{C}_{m} + { }_{n-1}\mathrm{C}_{m-1} = \frac{(n-1)!}{(n-m-1)!~m!}+\frac{(n-1)!}{(n-m)!~(m-1)!}\\
=(n-1)!~\frac{(n-m)+m}{(n-m)!~m!} = \frac{n!}{(n-m)!~m!} = { }_{n}\mathrm{C}_{m}$$

# 3. の証明
少し書き換えた
$${}_{m}\mathrm{C}_{m}+  {}_{m+1}\mathrm{C}_{m}+{}_{m+2}\mathrm{C}_{m}+...+{}_{m+n}\mathrm{C}_{m}= { }_{m+n+1}\mathrm{C}_{m+1}$$
を数学的帰納法で証明する。

i) n = 0のとき、
$${}_{m}\mathrm{C}_{m}= { }_{m+1}\mathrm{C}_{m+1}$$
より明らか。

ii) n = k が成り立つと仮定すると、n = k + 1のとき、
$${}_{m}\mathrm{C}_{m}+  {}_{m+1}\mathrm{C}_{m}+{}_{m+2}\mathrm{C}_{m}+...+{}_{m+k}\mathrm{C}_{m}+{}_{m+k+1}\mathrm{C}_{m} \\
= {}_{m+k+1}\mathrm{C}_{m+1} + {}_{m+k+1}\mathrm{C}_{m}$$
ここで、4. の等式を使うと、
$${}_{m+k+1}\mathrm{C}_{m+1} + {}_{m+k+1}\mathrm{C}_{m} = {}_{m+(k+1)+1}\mathrm{C}_{m+1}$$
よって、n = k + 1 でも成り立つ。

これより、$n \geq 0$ で 3. は示された。

グレた確率統計 ～負の二項分布再考～

前回は、畳み込みとそれを用いた再生性についてやりました。
二項分布には再生性がある、ということでしたね。
さて、それでは、幾何分布には再生性はないのでしょうか？

# 幾何分布と畳み込み

幾何分布に従う独立な確率変数２つの和を確率変数とする分布を求めてみましょう。

幾何分布は、

$$Pr(X=k)=pq^{k-1}$$

でしたね。

さて、独立な２つの確率変数、XとYを用意して、それらはそれぞれ

$$Pr(X=i)=pq^{i-1}$$

$$Pr(Y=j)=pq^{j-1}$$

とします。

そこで、新たな確率変数$Z=X+Y$を定義して、$Pr(Z=k)$ を求めてみようということです。

$$Pr(Z=k)= \sum_{l=1}^{k-1} Pr(X=l)Pr(Y=k-l) \\ = \sum_{l=1}^{k-1} pq^{l-1} \cdot pq^{k-l-1} = (k-1)p^2 q^{k-2}$$

どうやら、幾何分布には再生性は無いようです。

しかし、これは、負の二項分布のr=2に相当します！

どういうことかというと、新しい確率変数Zの意味するところが、

「事象がちょうど2回起こるのにかかる試行回数」だからです。

それでは、「事象がちょうど3回起こるのにかかる試行回数」は、というと、

$W = Z + G_1$ (ここで$G_1$は幾何分布に従う確率変数) なる $W$を定義すると、

$$Pr(W=k)= \sum_{l=1}^{k-1} Pr(W=l)Pr(G_1=k-l) \\ = \sum_{l=2}^{k-1} (l-1)p^2 q^{l-2} \cdot pq^{k-l-1} = p^3 q^{k-3} \frac{(k-2)(k-1)}{2}$$

となり、やはり、負の二項分布のr = 3 に相当します。

さあ、任意のrな負の二項分布に従う確率変数Xと、幾何分布に従う確率変数Yの和として定義した確率変数Zがr+1な負の二項分布に従うことを示しましょう。

負の二項分布は
$${}_{k-1}C_{r-1}~p^{r}(1-p)^{k-r}$$
でしたね。

$$Pr(Z=k)= \sum_{l=r}^{k-1} Pr(X=l)Pr(Y=k-l) \\ = \sum_{l=r}^{k-1} {}_{l-1}C_{r-1}~p^{r}q^{l-r} \cdot pq^{k-l-1} \\ =p^{r+1} q^{k-(r+1)} \sum_{l=r}^{k-1} {}_{l-1}C_{r-1} $$

ここで、
$$\sum_{l=r}^{k-1} {}_{l-1}C_{r-1}$$
は、$$ x = l -r$$と置くと、
$$\sum_{l=r}^{k-1} {}_{l-1}C_{r-1} = \sum_{x=0}^{k-r-1} {}_{x+r-1}C_{r-1}$$

組み合わせの公式に、
$${}_{m}\mathrm{C}_{m}+  {}_{m+1}\mathrm{C}_{m}+{}_{m+2}\mathrm{C}_{m}+...+{}_{n}\mathrm{C}_{m}= { }_{n+1}\mathrm{C}_{m+1}$$
があるので、これを用いると、結果的に、
$${ }_{k-1} \mathrm{C}_{r}~p^{r+1} q^{k-(r+1)}$$

が得られます。これはやはり、負の二項分布のr=r+1に相当しますね！

ということで、n個の幾何分布に従う確率変数の和はr=nの負の二項分布に従うわけです。

これは、モーメント母関数からも割と明らかで、
幾何分布に従う独立な確率変数、$X_1, X_2, ... , X_n$の和であるような確率変数Yのモーメント母関数は、
$$M_Y(t) = E[e^{t(X_1+X_2+X_3+...+X_n)}]=E[e^{tX_1}]E[e^{tX_2}]...=\\
\left(\frac{pe^t}{1-qe^t} \right )\left(\frac{pe^t}{1-qe^t} \right )\left(\frac{pe^t}{1-qe^t} \right )... =\left(\frac{pe^t}{1-qe^t} \right )^n$$
なわけですから、確かに幾何分布n個から成り立ってるなあと思えます。

ちなみに、二項分布に従う独立な確率変数、$X_1, X_2, ... , X_n$の和であるような確率変数Yのモーメント母関数は、
$$M_Y(t) = E[e^{t(X_1+X_2+X_3+...+X_n)}]=E[e^{tX_1}]E[e^{tX_2}]...\\
= (q+pe^t)^{n_1} (q+pe^t)^{n_2} (q+pe^t)^{n_3}... = (q+pe^t)^{n_1+n_2+n_3+...}$$
となるため、確かにモーメント母関数にも再生性が見られます。


このように、畳み込みを使うことで、より一般的な分布を導出することが可能なわけです。

# 少しの意味付け

上では、幾何分布を畳み込むことの動機を「再生性を確かめる」としましたが、

もう少し、応用的な側面で幾何分布を畳み込みたくなりましょう。

一旦、負の二項分布のことは忘れて（！）、幾何分布だけ知っている状態に戻りましょう！

すなわち、「初めて事象Hが起こるときの試行回数がk であるような確率」は知っているとします。

そして、「事象Hがちょうど2回起こるときの試行回数がkであるような確率」を求めてみましょう。

しかし、その前に、具体例を以って理解しておくことにします。

ex) 確率$p$で事象Hが起こるようなベルヌーイ試行がある。

このとき、事象Hが10回目でちょうど2回起こる確率はいくらか。

まず、この場合の幾何分布は、

$$Pr(X=k) = p \cdot q^{k-1}$$

となります。さて、10回目の時点でちょうど2回起こるというのは、

「10回目までのどこかで1回起こり、そして10回目でもう1回起こった」ということです。

「どこか」とはどこでもいいわけで、例えば、3回目で1回、１０回目でもう１回というような確率は、

$$ p q^2 \cdot p q^6$$

となります。ここで大事なのが、２回目が１０回目で生じたというのは、４回目を２度目の開始点とし、そこから数えて７回目にHが生じたということです。

「どこか」は他にもあって、

* １回目と１０回目（２回目から始めて９回目）

* ２回目と１０回目（３回目から始めて８回目）

* ３回目と１０回目（４回目から始めて７回目）

…

と全部で9つあります。これらをすべて足しあわせたものが答えになります。つまり、

$$ Pr (X=10) = \sum_{i=1}^9 p q^{i-1} \cdot p q^{(10-i)-1} $$

ここに畳み込みが現れているのが分かりますね。

この結果はr = 2 の負の二項分布で k=10 としたものと同じになるのは上で求めたとおりです。

これを踏まえると、単に10をkに置き換えれば、一般化できますね。

このように、ごく素朴に幾何分布を使って「k回目にちょうどr回起こる確率」を求めることができます。

グレた確率統計 ～畳み込みと再生性～

さて、前回は負の二項分布についてやりました。
負の二項定理なるものが出てきて、計算もやや煩雑で、うんざりしかけましたね。
今回は、畳み込みと、それを利用した再生性の問題について考えていきましょう。
結構重要ですが、本筋と関係ないといえば関係ありません。

# 畳み込み

畳み込みとは、簡単にいえば、

「$x + y = k$ の$x,y$を少しずつ変化させながら、足していく」という作業です。

なんじゃそりゃと思うかもしれませんが、確率の問題を解く際に無意識のうちに使っています。

ex1)

同様に確からしいサイコロ２つを投げる。その和が8となるような確率を求めよ。

ふつうに解いてみると、和が8となるような組み合わせは、

(2,6), (3,5), (4,4) の3通りで、サイコロは区別されるので、5/36 ですね。

少し丁寧な言い回しで解いてみましょう。

サイコロの出る目の確率変数をそれぞれ$X,Y$とします。

さて、今回求めるべき確率は、$Pr(X+Y = 8)$ですね。

XとYは独立ですから、これは、

$$\sum_{k=1}^6Pr(X=k)Pr(Y=8-k)$$

と書くことができます。はい、畳み込んでますよね！

一般的に、$Pr(X+Y = n)$を求めるとしても、やはり、

$$\sum_{k=1}^6Pr(X=k)Pr(Y=n-k)$$

とするだけです。これは、「２つのサイコロの目の和の確率分布」ですね。

そうすると、「２つのサイコロの目の和」という確率変数を導入したくなります。

これは、$Z = X+Y$ とおいてやればいいわけで、結局、

$$Pr(Z=n) = Pr(X+Y = n) = \sum_{k=1}^6Pr(X=k)Pr(Y=n-k)$$

ということです。

「複数の確率変数の和を新たな確率変数で表すとき」、畳み込みという操作がでてくるわけです。

# 二項分布の畳み込みと再生性

懐かしの二項分布に戻ってきましょう。
コインをn回投げて、表（H)の出る回数 X を調べていたとします。
そして、後日、また同じコインをm 回投げて、表の出る回数 Y を調べたとしましょう。
さて、この２つの実験を合わせた、表の出た回数 $i+j$ はどんな分布に従うでしょうか。

もちろん、$Z = X + Y$ と置きます。
また、それぞれの分布は、
$$Pr(X=i) = { }_n\mathrm{C}_i~p^iq^{n-i}\\
Pr(Y = j) = { }_m\mathrm{C}_j~p^jq^{m-j}$$
となります。さて、この２つの分布は独立ですから、
$$Pr(Z = n) = \sum_{k = 0}^n Pr(X=k)Pr(Y = n-k)$$
となります。
（総和の上限を n としたのは、別に $\infty$にしてもいいのですが、$k > n$ のときでは$Pr(Y=n-k）=0$ですから、結局足さなくてもいいわけです。）

ざくざく計算を進めてみましょう。
$$Pr(Z = k) = \sum_{l = 0}^k Pr(X=l)Pr(Y = k-l)\\
=\sum_{l = 0}^k { }_n\mathrm{C}_l~p^lq^{n-l} \cdot { }_m\mathrm{C}_{k-l}~p^{k-l}q^{m-(k-l)}\\
= p^k q^{n+m-k}\sum_{l = 0}^k { }_n\mathrm{C}_l \cdot { }_m\mathrm{C}_{k-l}$$
ここで、気をつけてほしいのが、 $k > n+m$ では確率は0であることです。
すなわち、上の計算は $0 \leq k \leq n+m$ に限っています。
さらに言うと、 $k > min(n,m)$ のとき、組み合わせCが未定義な項が生じます(${ }_nC_{n+1}$ とか)。
このような場合も確率は0となりますから、除外しています。
ですから、もう少し厳密に式を書くと、
$$ \mathrm{if}~k>n+m, Pr(Z=k)=0,\\
\mathrm{else}~Pr(Z=k) = p^k q^{n+m-k}\sum_{l = k-\mathrm{min}(k,m)}^{\mathrm{min}(k,n)} { }_n\mathrm{C}_l \cdot { }_m\mathrm{C}_{k-l} $$

さて、$\sum_{l = k-\mathrm{min}(k,m)}^{\mathrm{min}(k,n)} { }_n\mathrm{C}_l \cdot { }_m\mathrm{C}_{k-l}$ですが、
これは、$(n+m)$個の選択肢の中から$k$個選ぶときの場合の数は、選択肢を$n$個と$m$個に分け、$n$個のほうから$l$個、$m$個のほうから$k-l$個選ぶ、というのを$l$の取り得る範囲で行い、それぞれでの場合の数をすべて足すことによっても求められるということを意味しています。
きちんと証明したい人は、数学的帰納法で証明してください。読者への(ry

結局！
$$ \mathrm{if}~k < n+m, Pr(Z=k)=0,\\ \mathrm{else}~Pr(Z=k) = p^k q^{n+m-k}{ }_{n+m}\mathrm{C}_k$$
となるわけです。これは、二項分布の形ですね。

すなわち、独立な二項分布に従う２つの確率変数の和はやはり二項分布に従うということです。
このように独立な同分布に従う２つの確率変数の和が同じ分布に従うことを再生性があると言います。

グレた確率統計 ～負の二項分布～

さて、前回は二項分布と幾何分布についてやりました。
とはいえ、大本はベルヌーイ試行（0か1か）でしたね。
それぞれの意味合いは、
* 二項分布 : n回試行したときに、事象Hがk回起こる確率の分布
* 幾何分布 : 事象Hが初めて起こるのがk回目である確率の分布
でしたね。

ここから脇道に逸れるとすれば、多項分布もありますが、
ここでは伏線にもなりうる負の二項分布をやりたいと思います。

# 負の二項分布

負の二項分布は、幾何分布の拡張とも言えます。
すなわち、「事象Hがちょうどr回起こるのがk回目である確率」です。
これは、k-1 回目まででHが r-1 回起こり、k回目でHが起こる確率です。
少し考えれば容易に立式できて、
$$
Pr(X=k)= {}_{k-1}C_{r-1}~p^{r-1}(1-p)^{(k-1)-(r-1)}\cdot p \\
=  {}_{k-1}C_{r-1}~p^{r}(1-p)^{k-r}
$$
となります。
ただし、k < r のときは、Cが定義されないので、0とします。
計算上は、(k-1)回試行の二項分布で、k = r-1としたものに pをかけた形となっています。
もちろんのことながら、r = 1 とすると幾何分布になりますね。

形としては、二項分布と幾何分布の間の子…みたいな感じになりますね。

さて、期待値と分散を求めようと思うのですが、一旦立ち止まって考えてみましょう。
幾何分布、すなわち1回事象Hが起こる回数の期待値は$1/p$ でした。
単純に考えると、r回起こる回数の期待値は$r/p$のような気がしますね。
もちろん、負の二項分布の場合もモーメント母関数を使って、
平均と分散を出せるのですが、…天下り的に負の二項定理なる、
$$
\sum_{r=0}^\infty {}_{n+r-1}\mathrm{C}_{r}~x^r = (1-x)^{-n}
$$
を飲み込んでもらう必要があります。
これを飲み込んでもらうとして、モーメント母関数を求めると、
$$
M_X{(t)} = E[e^{tX}] = \sum_{k=r}^\infty {}_{k-1}C_{r-1}~e^{tk}p^{r}(1-p)^{k-r}
$$
ここで、$m = k-r$とおくと、${}_{s+t}\mathrm{C}_s = {}_{s+t}\mathrm{C}_t$、例の公式より、
$$
\sum_{k=r}^\infty {}_{k-1}C_{r-1}~e^{tk}p^{r}(1-p)^{k-r} \\
= \sum_{m=0}^\infty {}_{m+r-1}C_{r-1}~e^{t(m+r)}p^{r}(1-p)^{m}\\
= (pe^t)^r \sum_{m=0}^\infty {}_{m+r-1}C_{m}~(qe^t)^m\\
=(pe^t)^r\cdot(1-qe^t)^{-r}= \left( \frac{pe^t}{1-qe^t} \right )^r
$$
と求まります。なんと、これは幾何分布のモーメント母関数のr乗ですね！
幾何分布のモーメント母関数を$f_{(t)}$と表すことにすると、
$$
M_X(t)=(f{(t)})^r
$$
なわけですから、幾何分布のところで求めたモーメント母関数の計算結果を利用しつつ、

* 平均
$$
E[X] = \frac{\mathrm{d} }{\mathrm{d} t}M_X(t) \\
= r(f(t))^{r-1}\cdot f'(t)|_{t=0} = \frac{r}{p}
$$

* 分散
$$
E[X^2]=\frac{\mathrm{d^2} }{\mathrm{d} t^2}M_X(t) = \frac{\mathrm{d} }{\mathrm{d} t}\left( r f'(t)(f(t))^{r-1} \right)\\
= rf''(t)(f(t))^{r-1}+r(r-1)\left(f'(t) \right )^2 (f(t))^{r-2}|_{t=0} \\
=\frac{r(2-p)}{p^2}+\frac{r(r-1)}{p^2}
$$
$$
\therefore Var(X)=\frac{r(2-p)}{p^2}+\frac{r(r-1)}{p^2}-\left( \frac{r}{p} \right )^2 = \frac{r(1-p)}{p^2}
$$

となります。
予想通り（？）、平均も分散も幾何分布のr倍となっていますね！

※ 補足 ---- ここから ----
ある分布 f, gのモーメント母関数をそれぞれ $M_f(t),~M_g(t)$ とし、
$$
M_g(t)=(M_f(t))^n
$$
のようにn乗の関係にあるとき、
$$
E_g[X] = \frac{\mathrm{d} }{\mathrm{d} t}M_g(t) = n(M_f(t))^{n-1}\cdot {M_f}'(t)|_{t=0} = n{M_f}'(0) = E_f[X]
$$
$$
E_g[X^2]=\frac{\mathrm{d^2} }{\mathrm{d} t^2}M_g(t) = \frac{\mathrm{d} }{\mathrm{d} t}\left( n {M_f}'(t)({M_f}(t))^{n-1} \right)\\
= n{M_f}''(t)({M_f}(t))^{n-1}+n(n-1)\left({M_f}'(t) \right )^2 ({M_f}(t))^{n-2}|_{t=0}\\
= n{M_f}''(0) + n(n-1)\left({M_f}'(0) \right )^2
$$
$$
\therefore Var_g(X) = n{M_f}''(0) + n(n-1)\left({M_f}'(0) \right )^2 - \left(n{M_f}'(0)\right)^2\\
= n({M_f}''(0)-{M_f}'(0)) = n\cdot Var_f(X)
$$
と、平均、分散がそれぞれn倍になります。

---- ここまで ----

少し長くなったので、ここで一旦終わります。

グレた確率統計 ～二項分布と幾何分布～

（ここに、”確率とはなにか”というのがいるとは思うのですが、前提知識としておきます）
（期待値、分散、モーメント母関数の定義はwikipediaに任せます）

# 二項分布

「二項」と言うくらいですから、コインの裏表のように「２つの事象」のどちらかが生じるような試行を考えます。コインになぞらえて、２つの事象をそれぞれH(head)とT(tail)と呼ぶことにします。
各試行は独立とします。このような試行をベルヌーイ試行と言うそうです。
さて、ベルヌーイ試行をn回行ったとき、Hがk回生じる確率は、

$$
Pr(X = k) = {}_n \mathrm{C}_k p^k (1-p)^{n-k} = \frac{n!}{(n-k)! k!} p^k (1-p)^{n-k}
$$

となるのは高校出てたら分かるでしょう！
もちろん、これは $ k : 0 \to \infty $ まで足し合わせると、1になるのは読者への演習(ry

概形はおおよそこんな感じになります。

0.3の確率で当たるなら、10回繰り返せば3回くらい出るだろうというのは間違ってないのがわかります。



一応、期待値と分散を求めておきましょう。（以下、適宜 $ q = 1-p $と書きます）
そのために、モーメント母関数を求めると、

$$
M_X (t) = E[e^{tX}] = \sum_{k=1}^n {}_n\mathrm{C}_k~e^{tk}~p^k(1-p)^{n-k} \\ = \sum_{k=1}^n {}_n\mathrm{C}_k~(pe^{t})^k~q^{n-k} = (q +pe^t)^n
$$

なので、ここから芋づる式に期待値と分散は、

* 期待値
$$
E[X] = \frac{\mathrm{d} }{\mathrm{d} t}M_X(0) =npe^t(q+pe^t)^{n-1}|_{t=0} = np
$$

* 分散
$$
E[X^2] = \frac{\mathrm{d^2} }{\mathrm{d} t^2}M_X(0) = \frac{\mathrm{d} }{\mathrm{d} t} npe^t(q+pe^t)^{n-1} \\ = npe^t(q+pe^t)^{n-1}+n(n-1)p^2e^{2t}(q+pe^t)^{n-2} |_{t=0} \\ = np + n(n-1)p^2 = n^2p^2 - np^2 + np
$$
$$
\therefore Var[X] = E[X^2] - E[X]^2 \\ = n^2p^2 - np^2 + np - n^2p^2  = np(1-p)~~~(= npq)
$$

となりますね！

# 幾何分布

これもベルヌーイ試行を考えます。
先ほどの二項分布では、Hの生じた回数をカウントしていたわけですが、
今回は、Hが１回生じるまでにかかる試行回数をカウントしましょう。
試行の度にカウンターを押し、Hが出たときのカウンターの数値がk となっている確率は、
k-1 回Tが出て、その次にHが出ればいいわけですから、

$$
Pr(X=k)=p(1-p)^{k-1} ~~~ (= pq^{k-1})
$$

と簡単に求まります。これが $ k : 1 \to \infty $ まで足し合わせると１になるのは読者への(ry
概形はこんな感じになります。もちろんのことながら、右下がりの指数関数ですね。

これも期待値と分散を求めましょう。
もちろん！モーメント母関数を求めますと、

$$
M_X(t)=E[e^{tX}]= \sum_{k=1}^\infty e^{tk}p(1-p)^{k-1} \\
= \frac{p}{q}\sum_{k=1}^\infty (qe^t)^k = p\cdot \frac{e^t}{1-qe^t}
$$

となりますから、ここから期待値と分散は

* 期待値
$$
E[X] = \frac{\mathrm{d} }{\mathrm{d} t}M_X(0) = p\cdot\frac{e^t(1-qe^t)+qe^{2t}}{(1-qe^t)^2} = p\cdot \frac{e^t}{(1-qe^t)^2}|_{t=0}=\frac{1}{p}
$$

* 分散
$$
E[X^2] = \frac{\mathrm{d^2} }{\mathrm{d} t^2}M_X(0) = \frac{\mathrm{d} }{\mathrm{d} t} \left( p\cdot\frac{e^t}{(1-qe^t)^2}\right) \\
= p \cdot \frac{(1-qe^t)^2~e^t + 2qe^{2t}(1-qe^t)}{(1-qe^t)^4} \\
= p \cdot \frac{(1-qe^t)~e^t + 2qe^{2t}}{(1-qe^t)^3} |_{t=0} = p \cdot \frac {p+2q}{p^3}=\frac{2-p}{p^2}\\
$$
$$
\therefore
Var(X) = E[X^2]-E[X]^2 = \frac{2-p}{p^2} - \left( \frac{1}{p} \right )^2 = \frac{1-p}{p^2}
$$

となりますね。

さて、ここで幾何分布の無記憶性を示しておきましょう。
無記憶性というのは「今までの結果はこれからのことに影響しない」ということです。
数式で書くと、条件付き確率で

$$
\forall n,k \in \mathbb{N}~~Pr(X > n+k|X > n) = P(X > k)
$$

と書けます。これに代入するために、$Pr(X > k)$を求めておきましょう：

$$
Pr(X > k) = 1 - \sum_{i=1}^{k}p(1-p)^{i-1} \\
= 1 - p\cdot\frac{1-q^{k}}{1-q} = 1- (1-q^k) = q^k
$$

簡単ですね！てなわけで、無記憶性の数式に代入すると、

$$
Pr(X > n+k|X > n) = \frac{Pr(X > n+k \wedge  X>n)}{Pr(X > n)} \\
= \frac{Pr(X>n+k)}{Pr(X>n)} = \frac{q^{n+k}}{q^n} = q^k = Pr(X>k)
$$

となり、確かに「過去のことなんか、関係ない！」となります。

# どんなことに使えるか

これら２つの分布の大本はベルヌーイ試行です。つまり、「YES or NOな現象」です。
ですから、
* 繰り返し行われる
* 独立である
* 確率は一定である
* その結果がデジタル値(0か1か、YESかNOか、HかTか）
ならば、上２つの分布がなりたつと考えていいでしょう。

ex1)
商店街のイベントで電子くじ引きをすることに決まった。
電子くじ引きは当たりか外れかのどちらかが一定の確率で出るとする。
「どうせなら」ということで、当たりは一度出たらくじ引きは終了とし、
その代わり、景品を（しょぼくれた商店街のくせに）ハワイ旅行にすることに。
もちろんのことながら、すぐ当たられては商売上がったりである。
そこで、100人より多いところで当たりが出る確率が
i) 90%になる当たりの確率 p とそのときの期待値
ii) 95%になる当たりの確率 p とそのときの期待値
iii) 99% になる当たりの確率 p とそのときの期待値
を求めてほしい。

ex2)
商品のとある機械は、1ヶ月に1度のメンテナンスで一定の割合である部分の故障が見つかる。
うちの方針で、必ずこちらが修理代を負担することになっている。
i) 今月、50台の機械をメンテナンスしたところ、4台が故障していた。
これは、今までのデータを見る限りごく平均的な故障台数である。
機械１台が次の月に故障している確率 p を求めよ。

ii) こんなにしょっちゅう故障されては、修理代のせいで赤字になってしまう。
そこで、会議の末、「２年保証」にすることにした。
つまり、なんとか故障確率を下げて、２年まではもつようにすればいい。
そこで、２年経つまでは故障しない確率が90%になるような故障確率を求めてほしい。

iii) 改善が達成したとき、修理代が1回10万円かかるとして、
1台、1年あたりの経費削減の見込みを求めてほしい。

iv) さらに、現在の機械保有者数は50人であり、これから保有者数が変わらないと仮定する。
また、故障確率を0.001下げるのに20万円の費用がかかる。
経費削減の分でこれを賄おうとすると、大体何年かかるか。


～解答～
ex1)
i) p = 0.00105305, 期待値 : 950人
ii) p = 0.000512801, 期待値 : 1950人
iii) p = 0.000100498, 期待値 : 9950人

ex2)
i) p = 0.08
ii) $p_{\mathrm{new}} = 0.0044$
iii) 90720円
iv) 3.3年

グレた確率統計

統計学の本は大体、二項分布から始まり、ポアソン分布、幾何分布、指数分布、…と次に分布の一覧をよこしてきます。
ただ、分布をたくさん与えられる（押し付けられる？）からといって、確率分布への理解が深まるとはあまり思えません。
 たしかに、期待値や分散の計算問題としては優秀ではありますが…それだけのために、ねえ？

大まかに分けて、「確率と統計」は２つの道に進めると思います。
ひとつは、正統派な統計の方。推定とか検定に繋がる道です。こっちの道の主役は正規分布、t分布、χ2分布などなど…。実践的な推定や検定へとつながっていきます。
もうひとつは、「統計」らしからぬ（？）方で、指数分布、ガウス分布、ポアソン分布などが主役となります。 これはどんな道かというと、確率過程に繋がるであろう道です。待ち行列とかですね。 人によっては、統計の道よりもこの確率過程の道をガンガン進んだほうが楽しく、有用かもしれません。
（ツールとしての統計学においては、変に深く数学的基礎をやるよりは、さっさと推定・検定をエクセルとかRを使いながらざくざくやっていく方がいいかもしれないわけです。途中で詰まってしまうくらいなら！）

というわけで、おそらくどこかにそのようなテキストはあると思いますが、備忘録がてら、数回に分けてそれについて書いていこうかなと思うわけです。