ベルヌーイ分布を起点に確率分布間の関係をまとめてみた

# はじめに
　確率分布間の関係の理解を目的として、ベルヌーイ分布を起点に、様々な確率分布を生成する関数をPythonで実装し、まとめました。本記事は以下の方向けです。

* 基本的な確率分布は一通り勉強したことがある
* Pythonを触ったことがある
* 確率分布間の関係を理解したい

　実装したPythonの関数は以下の通りです。`bernolli`、`binom`などが関数名です。

![kankei1.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/571330/3ef28e31-3972-53c3-fda2-844657a600b5.png)

　例えば、一様分布は次の順にたどることで生成できます。
`ベルヌーイ分布 -> 幾何分布 -> 指数分布 -> ガンマ分布 -> ベータ分布 -> 一様分布`
　関数名や確率分布の定義は[scipy.stats](https://docs.scipy.org/doc/scipy/reference/stats.html)を踏襲しています。
　
　それでは、これから各関数を順に説明していきます。

# ベルヌーイ分布
　まず、出発点のベルヌーイ分布に従うサンプルを生成します。ベルヌーイ分布とは、「当たりが出る確率が$p$のくじを1回引いたときの当たり回数」に従う確率分布です。定義からも分かる通り、とりうる値は0か1です。確率関数は下式で与えられます。

```math
f(x;p) = p^x (1-p)^{1-x}\ \ (x=0,1)
```

<details><summary>　ベルヌーイ分布の平均は$p$、分散は$p(1-p)$です。</summary>

:::note info

[確率関数からの導出]

```math
\begin{align}
E(X) &= 0\times(1-p) + 1\times p = p　\\
V(X) &= E[(X-p)^2] \\ &= p^2 \times (1-p) + (1-p)^2 \times p \\&= p(1-p)
\end{align}
```
:::
</details>

　Python実装に進みます。あらかじめnumpyとscipy.statsをインポートしておきます。

```Python
import numpy as np
from scipy import stats
```

　確率が$p$のベルヌーイ分布に従うサンプルをsampleSize個発生させる関数は下の通りです。

```Python
def bernoulli(p, sampleSize):
    return stats.bernoulli.rvs(p, size=sampleSize)
```
試しにこの関数から20個のサンプルを生成してみます。


```Python
bernoulli(0.2, 20)
> array([1, 0, 0, 0, 1, 0, 0, 1, 0, 1, 0, 1, 0, 0, 0, 0, 0, 1, 0, 1])
```
設定どおり、約0.2の確率で1が出現しています。
　scipy.statsを使うのはこれが最後で、以後この`bernoulli関数`を元に確率分布を生成します。

# ベルヌーイ分布から二項分布
　二項分布とは、「当たりが出る確率が$p$のくじを$n$回引いた時の当たり回数」に従う確率分布です。確率関数は下式で与えられます。

```math
f(x;n,p) = \left(
\begin{matrix}
n\\
x
\end{matrix}
\right) p^x(1-p)^{n-x}\ \ (x=0,1,\cdots,n)
```

<details><summary>二項分布の平均は$np$、分散は$np(1-p)$です。</summary>

:::note info

[確率関数からの導出]

```math
\begin{align}
E(X) &= \sum_{x=0}^n x\left(\begin{matrix}
n\\
x
\end{matrix}\right)p^x(1-p)^{n-x} \\

&= np\sum_{x=1}^n \left(\begin{matrix}
n-1\\
x-1
\end{matrix}\right)p^{x-1}(1-p)^{n-x} \\

&= np\sum_{y=0}^{n-1} \left(\begin{matrix}
n-1\\
y
\end{matrix}\right)p^{y}(1-p)^{n-1-y} \\
&= np \\


E[X(X-1)] &=  \sum_{x=0}^n x(x-1) \left(\begin{matrix}
n\\
x
\end{matrix}\right)p^x(1-p)^{n-x} \\ 
&= 
n(n-1)p^2\sum_{x=2}^n \left(\begin{matrix}
n - 2\\
x - 2
\end{matrix}\right)p^{x-2}(1-p)^{n-x} \\

&= 
n(n-1)p^2\sum_{y=0}^{n-2} \left(\begin{matrix}
n - 2\\
y
\end{matrix}\right)p^{y}(1-p)^{n-2-y} \\
&=n(n-1)p^2 \\


V(X) &= E(X^2)-E(X)^2 \\ &= E[X(X-1)] + E(X) - E(X)^2 \\ &= n(n-1)p^2 + np - (np)^2 \\ &= np(1-p)
\end{align}

```

:::
</details>


　二項分布はベルヌーイ分布における「くじを引く回数」を$n$倍にしたものです。sampleSize個の「試行回数$n$、確率$p$の二項分布から生成されるサンプル」は、$n\times$sampleSize個の「確率$p$のベルヌーイ分布から生成されるサンプル」から作ることができます。実装は下の通りです。 


```Python
def binom(n, p, sampleSize):
    x = bernoulli(p, n * sampleSize)
    x = x.reshape([n, sampleSize])
    return np.sum(x, axis=0)
```


# 二項分布からポアソン分布
　ポアソン分布とは「単位時間あたり平均$\mu$個の当たりが出るくじにおいて、単位時間に出る当たりの数」に従う確率分布です。確率関数は下式で与えられます。

```math
f(x;\mu) = \frac{\mu^x}{x!}\exp(-\mu)\ \ (x=0,1,\cdots)
```

<details><summary>ポアソン分布の平均、分散はともに$\mu$です。</summary>

:::note
[確率関数からの導出]

```math
\begin{align}
E(X) &= \sum_{x=0}^\infty x\cdot\frac{\mu^x}{x!}\exp(-\mu) \\
&= \sum_{x=1}^\infty \cdot\frac{\mu\cdot\mu^{x-1}}{(x-1)!}\exp(-\mu) \\
&=\mu
\\

E[X(X-1)] &=  \sum_{x=0}^{\infty} x(x-1) \cdot \frac{\mu^x}{x!}\exp(-\mu) \\
&= \sum_{x=2}^\infty \cdot\frac{\mu^2\cdot\mu^{x-2}}{(x-2)!}\exp(-\mu) \\
&=\mu^2
\\

V(X) &= E(X^2)-E(X)^2 \\
&= E[X(X-1)] + E(X) - E(X)^2 \\
&= \mu^2 + \mu - \mu^2 \\
&= \mu
\end{align}

```
:::

</details>

下の通り、ポアソン分布は二項分布とよく似ています。

||二項分布|ポアソン分布|
|---|---|---|
|何を見る？|当たり回数|当たり回数|
|どういう条件で|試行回数一定|時間を一定|

試行回数$n$、確率$\frac{\mu}{n}$の二項分布の$n$を$n\to\infty$とすると、平均$\mu$のポアソン分布に収束します。実装は下の通りです。

```Python
def poisson(mu, sampleSize):
    n = 1000 # 十分大きい数
    p = mu / n
    return binom(n, p, sampleSize)
```

$n$は理論的には大きいければ大きいほど良いですが、計算量が増大します。
　ポアソン分布に従うサンプルは`binom関数`から生成しました。その`binom関数`は`bernoulli関数`が元になっています。つまり、ベルヌーイ分布からポアソン分布を生成したことになります。今後もこのような方法をとっていきます。

# ベルヌーイ分布から幾何分布
　幾何分布とは、「当たりが出る確率が$p$のくじを当たるまで引いたときの試行回数」に従う確率分布です。確率関数は下式で与えられます。

```math
f(x;p)= p(1-p)^{x-1}\ \ (x=1,2,\cdots)
```

<details><summary>幾何分布の平均は$\frac{1}{p}$、分散は$\frac{1-p}{p^2}$です。</summary>

:::note

[確率関数からの導出]

```math
E(X) = \sum_{x=1}^\infty x p(1-p)^{x-1} 
```
これを$S$とおく。
```math
\begin{align}
S &= \sum_{x=1}^\infty x p(1-p)^{x-1} \cdots(1)\\
  &= \sum_{x=0}^\infty (x+1) p(1-p)^x \\ 
  &= p+\sum_{x=1}^\infty (x+1) p(1-p)^x \cdots(2)\\ 
\end{align}
```

$(1-p)\times(1)$をすると、

```math
(1-p)S = \sum_{x=1}^\infty x p(1-p)^x \cdots(3)\\
```

$(2)-(3)$をすると

```math
pS=p+ \sum_{x=1}^\infty p(1-p)^x = 1
```

よって

```math
E(X)=S=\frac{1}{p}
```

次に、$E(X^2)$を求める

```math
E(X^2)=\sum_{x=1}^\infty x^2 p(1-p)^{x-1}
```

これを$T$とおくと

```math
\begin{align}
T&=\sum_{x=1}^\infty x^2 p(1-p)^{x-1} \cdots(4)\\
&=\sum_{x=0}^\infty (x+1)^2 p(1-p)^x \\
&=p+\sum_{x=1}^\infty (x+1)^2 p(1-p)^x \cdots(5)
\end{align}
```

$(1-p)\times(4)$をすると、

```math
(1-p)T=\sum_{x=1}^\infty x^2 p(1-p)^x \cdots(6)
```

$(5)-(6)$をすると、

```math
pT = p+\sum_{x=1}^\infty(2x+1)p(1-p)^x
```

$U=\sum_{x=1}^\infty(2x+1)p(1-p)^x$とおくと、

```math
\begin{align}
U&=\sum_{x=1}^\infty(2x+1)p(1-p)^x \cdots(7)\\
&=\sum_{x=0}^\infty(2x+3)p(1-p)^{x+1} \\
&=3p(1-p)+\sum_{x=1}^\infty(2x+3)p(1-p)^{x+1} \cdots(8) \\
\end{align}
```

$(1-p)\times(7)$をすると、

```math
(1-p)U=\sum_{x=1}^\infty(2x+1)p(1-p)^{x+1} \cdots(9)
```

$(8)-(9)$をすると、

```math
\begin{align}
pU&=3p(1-p)+\sum_{x=1}^\infty 2p(1-p)^{x+1}\\
&=3p(1-p)+2(1-p)^2\\
&=-p^2-p+2
\end{align}
```

以上より、

```math
\begin{align}
E(X^2)&=T=\frac{p+S}{p}=\frac{2-p}{p^2} \\
V(X)&=E(X^2)-E(X)^2= \frac{2-p}{p^2}-\frac{1}{p^2}=\frac{1-p}{p^2}
\end{align}
```
:::
</details>

確率$p$の幾何分布から生成されるサンプルは、確率$p$のベルヌーイ分布から作ることができます。実装は下の通りです。

```Python
def geom(p, sampleSize):
    x = []
    for _ in range(sampleSize):
        t = 1
        while bernoulli(1 - p, 1):
            t += 1
        x.append(t)
    return np.array(x)
```
```Python
def geom(p, sampleSize):
    n = 1000  
    x = []
    s = np.array([-1], dtype=int)
    while len(x) < sampleSize:
        x_ = bernoulli(p, n) # ベルヌーイ分布に従うサンプルを同時にnサンプル生成する。
        x_ = np.where(x_ == 1)[0]
        x_ = np.concatenate([s, x_])
        x.extend(x_[1:] - x_[:-1])
        s = np.array([x_[-1] - n])
    return np.array(x[:sampleSize])
```
　2つの関数は同じ処理をしていますが、前者のほうがわかりやすく、後者のほうが速いと思います。

# 幾何分布から負の二項分布
　負の二項分布とは、「当たりが出る確率が$p$のくじを$n$回当たるまで引いたときの、外れ回数」に従う確率分布です。確率関数は下式で与えられます。

```math
f(x;n,p) = \left(
\begin{matrix}n+x-1\\x\end{matrix}
\right) p^n(1-p)^{x}\ \ (x=0,1,\cdots,n)
```



<details><summary>負の二項分布の平均は$\frac{n(1－p)}{p}$、分散は$\frac{n(1-p)}{p^2}$です。</summary>

:::note info

[確率関数からの導出]

```math
\begin{align}
E(X) &= \sum_{x=0}^\infty x
\left(\begin{matrix}n+x-1\\x\end{matrix}\right)
p^n(1-p)^{x}\\
&= n\sum_{x=1}^\infty \left(\begin{matrix}n+x-1\\x-1\end{matrix}\right)p^n(1-p)^{x}\\
&= n\sum_{y=0}^\infty \left(\begin{matrix}n+y\\y\end{matrix}\right)p^n(1-p)^{y+1}\\
&= \frac{n(1-p)}{p}\sum_{y=0}^\infty \left(\begin{matrix}n+y\\y\end{matrix}\right)p^{n+1}(1-p)^{y}\\
&=\frac{n(1-p)}{p} \\

E[X(X-1)] &= \sum_{x=0}^\infty x(x-1)
\left(\begin{matrix}n+x-1\\x\end{matrix}\right)
p^n(1-p)^{x}\\
&=n(n+1)\sum_{x=2}^\infty 
\left(\begin{matrix}n+x-1\\x-2\end{matrix}\right)
p^n(1-p)^{x}\\
&=n(n+1)\frac{(1-p)^2}{p^2}\sum_{y=0}^\infty 
\left(\begin{matrix}n+y+1\\y\end{matrix}\right)
p^{n+2}(1-p)^{y}\\
&=\frac{n(n+1)(1-p)^2}{p^2}\\

V(X) &= E(X^2)-E(X)^2 \\
&= E[X(X-1)] + E(X) - E(X)^2 \\
&= \frac{n(n+1)(1-p)^2}{p^2} + \frac{n(1-p)}{p} - \frac{n^2(1-p)^2}{p^2}  \\
&= \frac{n(1-p)}{p^2}

\end{align}
```
:::
</details>

<details><summary>「負の二項分布」という名前の由来は、確率関数を変形すると二項分布に類似した形になるからです。

</summary>

:::note info

[確率関数の変形]

```math
\begin{align}
\left(\begin{matrix}n+x-1\\x\end{matrix}\right) p^n(1-p)^{x}&=\frac{(n+x-1)(n+x-2)\cdots(n+1)n}{x!}p^n(1-p)^x\\
&=\frac{-n(-n-1)\cdots(-n-x+1)}{x!}(-1)^xp^n(1-p)^x \\
&=\frac{-n(-n-1)\cdots(-n-x+1)}{x!}p^n(p-1)^x\\
&=\left(\begin{matrix}-n\\x\end{matrix}\right)p^n(p-1)^x
\end{align}
```
:::
</details>

　幾何分布との本質的な違いは、「当たりが何回出るまでくじを引くか」のみです。下の通り、負の二項分布は幾何分布から作ることができます。

```Python
def nbinom(n, p, sampleSize):
    x = geom(p, sampleSize * n) - 1
    x = x.reshape([sampleSize, n])
    return np.sum(x, axis=1)
```

# 幾何分布から指数分布
　指数分布とは、「単位時間あたり平均で$\frac{1}{\lambda}$個の当たりが出るくじにおいて、1回当たりが出るまでの時間」に従う確率分布です。確率密度関数は下式で与えられます。

```math
f(x,\lambda) = \lambda\exp(-\lambda x)
```



<details><summary>指数分布の平均は$\frac{1}{\lambda}$、分散は$\frac{1}{\lambda^2}$です。</summary>

:::note

[確率密度関数からの導出]

```math
\begin{align}
E(X)
&=\int_0^\infty \lambda x \exp(-\lambda x)dx \\
&= - \int_0^\infty x \{ \exp(-\lambda x)\}^{'} dx \\
&= \int_0^\infty  \exp(-\lambda x)dx \\
&= \frac{1}{\lambda} \\
E(X^2)
&=\int_0^\infty \lambda x^2 \exp(-\lambda x)dx \\
&= - \int_0^\infty x^2 \{ \exp(-\lambda x)\}^{'} dx \\
&= \int_0^\infty 2x \exp(-\lambda x) dx \\
&= \frac{2}{\lambda} E(X) \\
&= \frac{2}{\lambda^2} \\
V(X)
& = E(X^2) - E(X)^2\\
& = \frac{2}{\lambda^2} - \frac{1}{\lambda^2}\\
&= \frac{1}{\lambda^2}
\end{align}
```

:::

</details>


指数分布は次のように幾何分布とよく似ています。

||幾何分布|指数分布|
|---|---|---|
|何を見る？|1回当たりが出るまでの試行数|1回当たりが出るまでの時間|
|どういうくじ？|試行毎の当たり確率が一定|単位時間あたりの平均当たり回数が一定|

　時間や試行によらず、当たり確率が変わらない性質を**無記憶性**と呼びます。幾何分布も指数分布もこの性質を持っています。
　幾何分布において$p\to0$とすることで、指数分布を作ることができます。

```Python
def expon(lam, sampleSize):
    p = 0.01  # 十分小さい数
    x = geom(p, sampleSize)
    return x * p * lam
```

$p$は理論的には小さければ小さいほど良いですが、計算量が増大します。

# 指数分布からポアソン分布
　ポアソン分布は、二項分布だけでなく指数分布からも作ることができます。[指数分布とポアソン分布のいけない関係](https://www.slideshare.net/teramonagi/ss-11296227)のスライドが、非常に分かりやすいです。スライドにはRの実装が載っていますが、ここではPythonでの実装を示します。

```Python
def poisson_(mu, sampleSize):
    x = []
    while len(x) < sampleSize:
        t = 0
        n = -1
        while t <= 1:
            t += expon(1 / mu, 1)
            n += 1
        x.append(n)
    return np.array(x)
```

# 指数分布からガンマ分布
　ガンマ分布とは、「単位時間あたり平均で$\frac{1}{\lambda}$個当たりが出るくじにおいて、$\alpha$回当たりが出るまでの時間」に従う確率分布です。確率密度関数は下式で与えられます。

```math
f(x;\alpha,\lambda)=\frac{\lambda^\alpha}{\Gamma(\alpha)}x^{\alpha-1}\exp(-\lambda x)\ \ (x\ge0)
```

ただし、$\Gamma(\alpha)$は下式で与えられます。

```math
\Gamma(\alpha)=\int_0^\infty t^{\alpha-1}\exp(-t)\ dt
```


<details><summary>ガンマ分布の平均は$\frac{\alpha}{\lambda}$、分散は$\frac{\alpha}{\lambda^2}$です。</summary>

:::note

[確率密度関数からの導出]

```math
\begin{align}
E(X)
&=\int_0^\infty \frac{\lambda^\alpha}{\Gamma(\alpha)}x^{\alpha}\exp(-\lambda x)\ dx \\
&=\frac{\Gamma(\alpha+1)}{\lambda\Gamma(\alpha)}\int_0^\infty \frac{\lambda^{\alpha+1}}{\Gamma(\alpha+1)}x^{\alpha}\exp(-\lambda x)\ dx \\
&=\frac{\alpha}{\lambda} \\
E(X^2)
&=\int_0^\infty \frac{\lambda^\alpha}{\Gamma(\alpha)}x^{\alpha+1}\exp(-\lambda x)\ dx \\
&=\frac{\Gamma(\alpha+2)}{\lambda^2\Gamma(\alpha)}\int_0^\infty \frac{\lambda^{\alpha+2}}{\Gamma(\alpha+2)}x^{\alpha}\exp(-\lambda x)\ dx \\
&=\frac{\alpha(\alpha+1)}{\lambda^2}\\
V(X)&=E(X^2)-E(X)^2\\
&=\frac{\alpha}{\lambda^2}
\end{align}

```
:::
</details>

　指数分布との本質的な違いは、「当たりが何個出るまでの時間を見るか」のみです。ガンマ分布は指数分布から作ることができます。

```Python
def gamma(alpha, lam, sampleSize):
    x = expon(lam, sampleSize * alpha)
    x = x.reshape([sampleSize, alpha])
    return np.sum(x, axis=1)
```
尚、上記`gamma関数`の引数$\alpha$は自然数である必要があります。一般のガンマ分布の$\alpha$は正の数であればなんでもOKなのですが、自然数以外の$\alpha$だとベルヌーイ分布から作れないですし、説明も難しいので、ここでは自然数に限定しました。

# 負の二項分布からガンマ分布
　ガンマ分布は負の二項分布からも作ることができます。次のようにガンマ分布と負の二項分布はよく似ています。

||負の二項分布|ガンマ分布|
|---|---|---|
|何を見る？|$n$回当たりが出るまでの試行数|$\alpha$回当たりが出るまでの時間|
|どういうくじ？|試行毎の当たり確率が一定|単位時間あたりの平均当たり回数が一定|

　負の二項分布からガンマ分布を作る方法は、幾何分布から指数分布を作る方法と全く同じです。

```Python
def gamma_(alpha, lam, sampleSize):
    p = 0.01  # 十分小さい数
    x = nbinom(alpha, p, sampleSize)
    return x * p * lam
```

ここでも、引数$\alpha$は自然数である必要があります。

# 二項分布から正規分布
　どのような分布も無限に足し合わせると正規分布になります（中心極限定理）。私が初めて作った記事「[確率分布の漸近的性質をPythonで確かめる](https://qiita.com/yutera12/items/e9dfe9f6aafe96f3c700)」に例があります。
　二項分布はベルヌーイ分布の$n$個の足し合わせでした。$n\to\infty$とすることで二項分布は正規分布に従います。正規分布の確率密度関数は下式で与えられます。

```math
f(x;\mu,\sigma^2)=\frac{1}{\sqrt{2\pi\sigma^2}}\exp\biggl\{-\frac{(x-\mu)^2}{2\sigma^2}\biggr\}
```

<details><summary>正規分布の平均は$\mu$、分散は$\sigma^2$です。</summary>

:::note

[確率密度関数からの導出]

```math
\begin{align}
E(X)&=\int_{-\infty}^{\infty}\frac{x}{\sqrt{2\pi\sigma^2}}\exp\biggl\{-\frac{(x-\mu)^2}{2\sigma^2}\biggr\}\ dx \\
&=\int_{-\infty}^{\infty}\frac{y+\mu}{\sqrt{2\pi\sigma^2}}\exp\biggl\{-\frac{y^2}{2\sigma^2}\biggr\}\ dy\\
&=\int_{-\infty}^{\infty}\frac{y}{\sqrt{2\pi\sigma^2}}\exp\biggl\{-\frac{y^2}{2\sigma^2}\biggr\}\ dy + \mu\\
&=\mu \\
V(X)&=\int_{-\infty}^{\infty}\frac{(x-\mu)^2}{\sqrt{2\pi\sigma^2}}\exp\biggl\{-\frac{(x-\mu)^2}{2\sigma^2}\biggr\}\ dx \\
&=\int_{-\infty}^{\infty}\frac{y^2}{\sqrt{2\pi\sigma^2}}\exp\biggl\{-\frac{y^2}{2\sigma^2}\biggr\}\ dx \\
&=-\sigma^2 \int_{-\infty}^{\infty}\frac{y}{\sqrt{2\pi\sigma^2}}\Biggl[\exp\biggl\{-\frac{y^2}{2\sigma^2}\biggr\}\Biggr]^{'}\ dx \\
&=\sigma^2 \int_{-\infty}^{\infty}\frac{1}{\sqrt{2\pi\sigma^2}}\exp\biggl\{-\frac{y^2}{2\sigma^2}\biggr\}\ dx \\
&=\sigma^2
\end{align}
```
:::
</details>

平均$\mu$、標準偏差$\sigma$の正規分布に従うサンプルは、次の方法で生成することができます。

```Python
def norm(mu, sigma, sampleSize):
    n = 1000 # 十分大きい数
    p = 0.5
    x = binom(n, p, sampleSize)
    sd = (n * p * (1 - p)) ** 0.5
    x = (x - n * p) / sd * sigma + mu
    return x
```

# 正規分布からカイ二乗分布
　くじのような説明はできませんが、仮説検定でよく出てくるのが$\chi^2$分布です。$X_1,X_2,\cdots,X_{k}$が、独立に標準正規分布（平均0、分散1の正規分布）に従うとき、$X_1^2+X_2^2+\cdots+X_{k}^2$を自由度$k$の$\chi^2$分布と言います。
 カイ二乗分布の確率密度関数は下式で与えられます。

```math
f(x;k)=\frac{1}{2^{k/2}\Gamma(k/2)}x^{k/2-1}e^{-x/2}
```

<details><summary>カイ二乗分布の平均は$k$、分散は$2k$です。</summary>

:::note

[確率密度関数からの導出]

```math
\begin{align}
E(X)&=\int_{0}^{\infty}\frac{x}{2^{k/2}\Gamma(k/2)}x^{k/2 \ -1}e^{-x/2} \ dx\\
&=\frac{2\Gamma(k/2 \ +1)}{\Gamma(k/2)}\int_{0}^{\infty}\frac{1}{2^{k/2 \ +1}\Gamma(k/2 \ +1)}x^{k/2}e^{-x/2} \ dx\\
&=2\times\frac{k}{2}=k\\
E(X^2)&=\int_{0}^{\infty}\frac{x^2}{2^{k/2}\Gamma(k/2)}x^{k/2 \ -1}e^{-x/2} \ dx\\
&=\frac{2^2 \ \Gamma(k/2 \ +2)}{\Gamma(k/2)}\int_{0}^{\infty}\frac{1}{2^{k/2 \ +2}\Gamma(k/2 \ +2)}x^{k/2 \ +1}e^{-x/2} \ dx\\
&=2^2\times\Bigl(\frac{k}{2}+1\Bigr)\frac{k}{2}=k(k+2)\\
V(X)&=E(X^2)-E(X)^2=2k

\end{align}
```
:::
</details>
 
カイ二乗分布に従うサンプルは、次の方法で生成することができます。

```Python
def chi2(k, sampleSize):
    x = norm(0, 1, sampleSize * k) ** 2
    x = x.reshape([sampleSize, k])
    return np.sum(x, axis=1)
```

# カイ二乗分布からガンマ分布

　これも、くじのような説明できませんが、自由度1の$\chi^2$分布は、$\alpha=1/2$、$\lambda=1/2$のガンマ分布に一致します。ガンマ分布の**再生成**という性質より、$\alpha$が$\frac{1}{2}$の$n$倍（$n$は自然数）の場合のガンマ分布を生成できます。再生性については、最後に説明します。

```Python
def gamma__(alpha, lam, sampleSize):
    df = int(np.round(alpha * 2))
    x = chi2(1, sampleSize * df)
    x = x.reshape([sampleSize, df])
    x = x * lam / 2
    return np.sum(x, axis=1)
```

この関数の引数$\alpha$は$\frac{1}{2}$の$n$倍（$n$は自然数）である必要があります。

# ガンマ分布からベータ分布

　これもくじのような説明ができませんが、「パラメータが$(\alpha,\lambda)$のガンマ分布」と「パラメータが$(\beta,\lambda)$のガンマ分布」をもとに、パラメータが$(\alpha,\beta)$のベータ分布を生成することができます。
　ベータ分布の密度関数は下式で与えられます。

```math
f(x;\alpha,\beta)=\frac{x^{\alpha-1}(1-x)^{\beta-1}}{B(\alpha, \beta)}\ \ (0\le x\le 1)
```

ただし、$B(\alpha, \beta)$は下式で与えられます。

```math
B(\alpha, \beta)=\int_0^1 t^{\alpha-1}(1-t)^{\beta-1} dt
```


<details><summary>ベータ分布の平均は$\frac{\alpha}{\alpha+\beta}$、分散は$\frac{\alpha\beta}{(\alpha+\beta)^2(\alpha+\beta+1)}$です。</summary>

:::note

[確率密度関数からの導出]

```math
\begin{align}
E(X)
&=\int_0^1 \frac{x^{\alpha}(1-x)^{\beta-1}}{B(\alpha, \beta)} dx\\
&=\int_0^1 \frac{x^{(\alpha + 1)-1}(1-x)^{\beta-1}}{B(\alpha + 1, \beta)} dx \ \ \frac{{B(\alpha + 1, \beta)}}{{B(\alpha , \beta)}}\\
&=\frac{{B(\alpha + 1, \beta)}}{{B(\alpha , \beta)}}\\
&=\frac{\alpha}{\alpha+\beta}\\

E(X^2)
&=\int_0^1 \frac{x^{\alpha+1}(1-x)^{\beta-1}}{B(\alpha, \beta)} dx\\
&=\int_0^1 \frac{x^{(\alpha + 2)-1}(1-x)^{\beta-1}}{B(\alpha + 2, \beta)} dx \ \ \frac{{B(\alpha + 2, \beta)}}{{B(\alpha , \beta)}}\\
&=\frac{{B(\alpha + 2, \beta)}}{{B(\alpha , \beta)}}\\
&=\frac{\alpha+1}{\alpha+\beta+1}\frac{{B(\alpha + 1, \beta)}}{{B(\alpha , \beta)}}\\
&=\frac{\alpha(\alpha+1)}{(\alpha+\beta)(\alpha+\beta+1)}\\

V(X)&=E(X^2)-E(X)^2\\
&=\frac{\alpha\beta}{(\alpha+\beta)^2(\alpha+\beta+1)}
\end{align}

```
:::
</details>
取りうる範囲が0～1であるため、ベイズ統計学において、確率を表すパラメータの事前分布としてよく用いられます。


```Python
def beta(alpha, beta, sampleSize):
    x1 = gamma(alpha, 1, sampleSize)
    x2 = gamma(beta, 1, sampleSize)
    return x1 / (x1 + x2)
```
# ベータ分布から一様分布
区間$[a,b]$の一様分布の密度関数は下式で与えられます。

```math
f(x;a,b)=\frac{1}{b-a}\ \ (b\le x\le a)
```
 <details><summary>　区間 $[a,b]$の一様分布の平均は$\frac{a+b}{2}$、分散は$\frac{(b-a)^2}{12}$です。</summary>

:::note info

[確率関数からの導出]

```math
\begin{align}
E(X) &= \int_a^b \frac{x}{b-a} dx = \frac{a+b}{2}\\
E(X^2) &= \int_a^b \frac{x^2}{b-a} dx = \frac{a^2+ab+b^2}{3}\\
V(X) &= E[X^2] - E[X]^2 = \frac{(b-a)^2}{12}
\end{align}
```
:::
</details>
 
 ベータ分布の密度関数から明らかなとおり、$\alpha,\beta＝1$の場合、ベータ分布は区間[0,1]の一様分布と一致します。したがって、一様分布に従うサンプルは、ベータ分布から作ることができます。


```Python
def uniform(sampleSize):
    return beta(1, 1, sampleSize)
```

# まとめ
![kankei2.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/571330/29577e08-257d-fda3-85cb-12c6f300ffc3.png)



　互いの関係が一目瞭然です。

　ここで、**再生性**について説明します。再生性とは分布の足し算を、パラメータの足し算に置き換えることができる性質のことです。例えば、

1. 試行回数が$n$、確率$p$の二項分布をX
1. 試行回数が$m$、確率$p$の二項分布をY

とすると、$X+Y$は、試行回数が$(n+m)$、確率$p$の二項分布に従います。1.は、$n$個のベルヌーイ分布の和、2.は$m$個のベルヌーイ分布の和に等しいので、$X+Y$は$(n+m)$個のベルヌーイ分布の和、すなわち、試行回数が$(n+m)$の二項分布と一致するからです。
　同じ理論で、上図の`複数化`の矢印の先の分布には、必ず再生性の性質があります（二項分布、負の二項分布、ガンマ分布、$\chi^2$分布）。そして、二項分布の極限である、正規分布とポアソン分布にも当然、再生性があります。

# 最後にPython実験
`ベルヌーイ分布 -> 幾何分布 -> 指数分布 -> ガンマ分布 -> ベータ分布 -> 一様分布 -> ベルヌーイ分布`を試してみます。

　あらかじめプロットのための関数を定義しておきます。

```Python
import matplotlib.pyplot as plt
import seaborn as sns
sns.set(style='whitegrid')

def plot_descreteDitribution(sample, truePf, xmin, xmax):
    sampleSize = len(sample)
    x = np.arange(xmin, xmax + 1)
    pf = np.array([np.sum(sample == i) for i in x]) / sampleSize  # 実験で求めた分布    
    
    # 実験で求めた分布を青色でプロット
    plt.plot(x - 0.1, pf, 'bo', ms=8)
    plt.vlines(x - 0.1, 0, pf, colors='b', lw=5, alpha=0.5)

    # 真の分布を赤色でプロット
    plt.plot(x + 0.1, truePf(x), 'ro', ms=8)
    plt.vlines(x + 0.1, 0, truePf(x), colors='r', lw=5, alpha=0.5)

def plot_continuousDitribution(sample, truePf, xmin, xmax):
    sampleSize = len(sample)
    x = np.linspace(xmin, xmax, 100)

    # 実験で求めた分布を青色でプロット
    th = np.linspace(xmin, xmax, 30)
    hi = np.array([np.sum(np.logical_and(th[i] < sample, sample <= th[i + 1])) for i in range(30 - 1)]) / (sampleSize * (th[1] - th[0]))
    plt.bar((th[:-1] + th[1:]) / 2, hi, width=(th[1] - th[0]))

    # 真の分布を赤色でプロット
    plt.plot(x, truePf(x), 'r', linewidth=4)

    plt.xlim(xmin, xmax)
```
#### 幾何分布(from ベルヌーイ分布)
 
```Python
p = 0.2
xmin = 1
xmax = 20

sample = geom(p, sampleSize)
truePf = lambda x: stats.geom.pmf(x, p)
plot_descreteDitribution(sample, truePf, xmin, xmax)
```

- 結果

![幾何分布.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/571330/4c53ba62-f5a8-d89d-b19c-ba5db1c3e0c2.png)

赤色がscipy.statsで作った分布、青色が今回作成した関数から作った分布です。尚、sampleSizeは10000としました。以後も同様です。

#### 指数分布(from ベルヌーイ分布 -> 幾何分布)

```Python
lam = 2
xmin = 0
xmax = 10

sample = expon(lam, sampleSize)
truePf = lambda x: stats.expon.pdf(x, scale=lam)
plot_continuousDitribution(sample, truePf, xmin, xmax)
```

- 結果

![指数分布.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/571330/52a6b2ba-195e-6fbc-1386-a3cd2fd404cf.png)

#### ガンマ分布(from ベルヌーイ分布 -> 幾何分布 -> 指数分布)

```Python
alpha = 3
lam = 2
xmin = 0
xmax = 20

sample = gamma(alpha, lam, sampleSize)
truePf = lambda x: stats.gamma.pdf(x, alpha, scale=lam)
plot_continuousDitribution(sample, truePf, xmin, xmax)
```

- 結果

![ガンマ分布.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/571330/b92b21e0-331b-f9fd-550a-78f7c65ba810.png)


#### ベータ分布(from ベルヌーイ分布 -> ... -> ガンマ分布)

```Python
alpha = 3
bet = 2
xmin = 0
xmax = 1

sample = beta(alpha, bet, sampleSize)
truePf = lambda x: stats.beta.pdf(x, alpha, bet)
plot_continuousDitribution(sample, truePf, xmin, xmax)
```

- 結果

![ベータ分布.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/571330/fe0da76e-17a9-8c64-d179-b5fc52e6ecf7.png)

#### 一様分布(from ベルヌーイ分布 -> ... -> ベータ分布)

```Python
xmin = 0
xmax = 1

sample = uniform(sampleSize)
truePf = lambda x: stats.uniform.pdf(x)
plot_continuousDitribution(sample, truePf, xmin, xmax)
```

- 結果

![一様分布.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/571330/c61822e0-fa64-8c87-a9db-f59e698b0bfa.png)

#### ベルヌーイ分布(from ベルヌーイ分布 -> ... -> 一様分布)！？
``` Python
p = 0.2
xmin = -2
xmax = 4

def bernoulli_(p, sampleSize):
    x = uniform(sampleSize)    
    return (x < p).astype(int)

sample = bernoulli_(p, sampleSize)
truePf = lambda x: stats.bernoulli.pmf(x, p)
plot_descreteDitribution(sample, truePf, xmin, xmax)
```

- 結果

![ベルヌーイ分布.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/571330/44be09ae-2606-ee70-577c-818314e3301e.png)

上手くいったことがわかります。

# 参考文献
- [色々な確率分布とその関係](https://www.slideshare.net/hirokiiida165/ss-78477986)
- [代表的な確率分布の特徴まとめ](https://qiita.com/qiita_kuru/items/d9782185652351c78aac)
- [scipy.stats](https://docs.scipy.org/doc/scipy/reference/stats.html)
- [確率分布の漸近的性質をPythonで確かめる](https://qiita.com/yutera12/items/e9dfe9f6aafe96f3c700)
- [指数分布とポアソン分布のいけない関係](https://www.slideshare.net/teramonagi/ss-11296227)
- [統計検定1級に合格する方法](https://qiita.com/drken/items/089b8443305df047b44e)