1. Pythonで学ぶ統計学　2-1. 確率分布［離散型変数］

*   **離散型の確率変数**は、サイコロの目のように飛び飛びの値をとる変数のことで、例えば｢1｣の次は｢2｣、｢2｣の次は｢3｣というように、その間に 1.1, 1.2, 1.3, ･･･, 1.8, 1.9 などといった連続的な数値は存在しません。
*   主な**離散型の確率分布**について、scipy.statsの`pmf`（probability mass function : 確率質量関数)や`rvs`（random variates : 確率変数）を使ってその特徴を見ていきます。

```python
# 数値計算ライブラリのインポート
import numpy as np
import scipy as sp
import pandas as pd
from pandas import Series, DataFrame
# 可視化ライブラリのインポート
import matplotlib.pyplot as plt
import matplotlib as mpl
import seaborn as sns
%matplotlib inline
# matplotlibの日本語表示モジュール
!pip install japanize-matplotlib
import japanize_matplotlib
```

##⑴ ベルヌーイ分布（Bernoulli distribution）
*   2種類しかない事象のどちらかを発生させることを**ベルヌーイ試行（Bernoulli trial）**といいます。
*   **ベルヌーイ分布**とは1回のベルヌーイ試行において、それぞれの事象が生じる確率分布のことです。
*   例えば、コインを8回投げて表が出たら0、裏が出たら1として結果を次のように仮定します。

```python
x = np.array([0,0,1,1,0,1,0,0])

# 確率分布を計算
p = len(x[x==1]) / len(x)
pmf_bernoulli = sp.stats.bernoulli.pmf(x, p)

# 可視化
plt.vlines(x, 0, pmf_bernoulli, 
           colors='blue', lw=50)
plt.xticks([0,1])
plt.xlim([0 - 0.5, 1 + 0.5])
plt.grid(True)
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/502109/1e0dab5c-c434-f830-0f21-0c1d35f06343.png)

|2種類の事象|確率|
|:---:|:---:|
|**0**|0.625|
|**1**|0.375|

##⑵ 二項分布（Binomial distribution）
*   互いに独立したベルヌーイ試行をn回くり返したとき、ある事象が何回起こるかという確率分布のことを**二項分布**といいます。
*   例えば、表が出る確率pが50％のコインを5回投げて、そのうち2回が表になる確率を`binom.pmf`を使って求めます。

```python
sp.stats.binom.pmf(n=5, p=0.5, k=2)
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/502109/5c2fd459-c4b0-f6af-beb9-4193330a1ca3.png)

*   表が出る確率pが20％のコインを10回投げて表が出た回数を数える、という試行を10000回くり返したと仮定し、`binom.rvs`を使って二項分布に従う疑似乱数を生成します。
*   さらに`binom.pmf`を使って、表が出る確率pが20％のコインを10回投げたときに表が出る回数の確率分布を計算し、疑似乱数のヒストグラムと比較してみます。

```python
# 疑似乱数を生成
np.random.seed(1)
rvs_binom = sp.stats.binom.rvs(n=10, p=0.2, size=10000)

# 確率分布を取得
m = np.arange(0, 10+1, 1)
pmf_binom = sp.stats.binom.pmf(n=10, p=0.2, k=m)

# 可視化
sns.distplot(rvs_binom, bins=m, 
             kde=False, norm_hist=True, label='rvs')
plt.plot(m, pmf_binom, label='pmf')
plt.xticks(m)
plt.legend()
plt.grid()
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/502109/be4f0a03-6354-e272-0240-9ce365c5e945.png)

|表が出る回数|確率|
|:---:|:---:|
|**0**|0.107374182|
|**1**|0.268435456|
|**2**|0.301989888|
|**3**|0.201326592|
|**4**|0.088080384|
|**5**|0.026424115|
|**6**|0.005505024|
|**7**|0.000786432|
|**8**|0.000073728|
|**9**|0.000004096|
|**10**|0.000000102|

##⑶ ポアソン分布（Poisson distribution）
*   ごく稀に起こる事象の確率で、例えば単位面積当たりの雨粒の数とか、ある交差点で一年間に発生する事故の回数などが**ポアソン分布**に従います。
*   つまり、一定の時間や面積に対して一定の割合で発生する確率分布で、サンプル数nが十分に大きく確率pが非常に小さい場合です。
*   例えば、ある期間に平均5回起こる現象が2回だけ起こる確率を`poisson.pmf`を使って求めます。

```python
sp.stats.poisson.pmf(k=2, mu=5)
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/502109/c1cd6cfc-556e-463d-aa8f-466aede5e3a2.png)

*   ポアソン分布のパラメータはある事象の**平均発生回数mu**で、これは**強度**とか**λ（ラムダ）**とも呼ばれます。
*   ある事象が起こる確率pを20％として、起こった回数を数えるという試行を10000回くり返したと仮定し、`poisson.rvs`を使ってポアソン分布に従う疑似乱数を生成します。
*   さらに`poisson.pmf`を使って、発生確率pが20％のときの確率分布を計算し、疑似乱数のヒストグラムと比較してみます。

```python
# 疑似乱数を生成
np.random.seed(1)
rvs_poisson = sp.stats.poisson.rvs(mu=2, size=10000)

# 確率分布を取得
m = np.arange(0, 10+1, 1)
pmf_poisson = sp.stats.poisson.pmf(mu=2, k=m)

# 可視化
sns.distplot(rvs_poisson, bins=m, 
             kde=False, norm_hist=True, label='rvs')
plt.plot(m, pmf_poisson, label='pmf')
plt.xticks(m)
plt.legend()
plt.grid()
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/502109/b1c032c7-2682-e780-14f5-629fbb83addd.png)

|発生回数|確率|
|:---:|:---:|
|**0**|0.135335283|
|**1**|0.270670566|
|**2**|0.270670566|
|**3**|0.180447044|
|**4**|0.090223522|
|**5**|0.036089409|
|**6**|0.012029803|
|**7**|0.003437087|
|**8**|0.000859272|
|**9**|0.000190949|
|**10**|0.000038190|

*   ここで、**ポアソン分布**と**二項分布**の関係を考察します。
*   試行回数nが十分に大きく確率pが非常に小さいときの**二項分布**の確率分布を計算し、先の**ポアソン分布**の例と比較してみます。

```python
# パラメータを指定
n = 100000000
p = 0.00000002

# 二項分布の確率分布を計算
num = np.arange(0, 10+1, 1)
pmf_binom_2 = sp.stats.binom.pmf(n=n, p=p, k=num)

# 可視化
plt.plot(m, pmf_poisson, 
         color='lightgray', lw=10, label='poisson')
plt.plot(m, pmf_binom_2, 
         color='black', linestyle='dotted', label='binomial')
plt.xticks(num)
plt.legend()
plt.grid()
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/502109/aaa31b7e-e7fe-964e-34fc-725762fe8a6e.png)

*   この場合、**ポアソン分布**と**二項分布**はほぼ一致し、両者が近い関係にあることがわかります。
*   **ポアソン分布**は、**二項分布**の発生確率pが非常に小さく試行回数nが十分に大きいときの状況を近似的に表します。

## ⑷ 幾何分布（geometric distribution）
*   成功確率がpの独立したベルヌーイ試行をくり返すとき、**初めて成功するまでの試行回数kが従う確率分布**を**幾何分布**といいます。
*   例えば、サイコロを１回だけ投げて「1」が出る確率をscipy.statsの`geom.pmf`を使って取得します。

```python
%precision 3
sp.stats.geom.pmf(k=1, p=1/6)
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/502109/3b3720ff-b4e5-dfff-65c3-f5f89afc572a.png)

*   サイコロの目は全部で6つ、すべてが同じ確率とすれば 1/6 で 0.167 になるのは当然ですが、
*   次々と投げて、2回目に初めて｢1｣が出る確率、3回目に初めて｢1｣が出る確率、･･････、10回目に初めて｢1｣が出る確率を求めます。

```python
# 試行回数を指定
num = np.arange(1, 11, 1)

# 確率分布を計算
prob = []
for i in num:
    value = sp.stats.geom.pmf(k=i, p=1/6)
    prob.append(value)

# 可視化
plt.bar(num, prob)
plt.xticks(num)
plt.xlabel('初めて1が出るまでの回数')
plt.ylabel('確率')
plt.show()
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/502109/17045120-223f-9834-15b1-5ba966794781.png)

*   k回投げて、k-1回目までは失敗し、k回目に初めて成功する確率ということで、この場合の確率分布は次の通りです。

|試行回数|確率|計算式|
|:---:|:---:|---|
|**1**|0.167|⅙|
|**2**|0.139|⅚・⅙|
|**3**|0.116|⅚・⅚・⅙|
|**4**|0.096|⅚・⅚・⅚・⅙|
|**5**|0.080|⅚・⅚・⅚・⅚・⅙|
|**6**|0.067|⅚・⅚・⅚・⅚・⅚・⅙|
|**7**|0.056|⅚・⅚・⅚・⅚・⅚・⅚・⅙|
|**8**|0.047|⅚・⅚・⅚・⅚・⅚・⅚・⅚・⅙|
|**9**|0.039|⅚・⅚・⅚・⅚・⅚・⅚・⅚・⅚・⅙|
|**10**|0.032|⅚・⅚・⅚・⅚・⅚・⅚・⅚・⅚・⅚・⅙|

##⑸ 離散一様分布（uniform distribution）
*   一様分布には離散型と連続型があり、確率変数が離散型の一様分布を**離散一様分布**といいます。
*   すべての事象の起こる確率が等しい分布で、例えばサイコロは1から6までの目が出る確率がすべて等しいので**離散一様分布**に従います。

```python
# 全事象を指定
num = np.arange(1, 7, 1)

# 確率分布を計算
prob = []
for i in num:
    value = 1 / len(num)
    prob.append(value)

# 可視化
plt.bar(num, prob)
plt.xticks(num)
plt.xlabel('サイコロの目')
plt.ylabel('確率')
plt.show()
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/502109/835d43e3-c05d-d1a5-6b2f-d8b590b06c1e.png)

|サイコロの目|確率|
|:---:|:---:|
|**1**|0.167|
|**2**|0.167|
|**3**|0.167|
|**4**|0.167|
|**5**|0.167|
|**6**|0.167|

##⑹ 超幾何分布（hypergeometric distribution）
*   例えば、全部で20個のくじがあり、そのうち7個が当たりであるとします。20個からランダムに12個を選んだとき、当たりが何個出るか。
*   この「引き当てる個数」を確率変数とする分布を**超幾何分布**といいます。
*   パラメータは3つで、全体の個数M、当たりの個数n、選択する個数Nになります。

```python
# パラメータを指定
M = 20 #全体個数
n = 7  #当たり個数
N = 12 #選択個数

# 確率変数を作成
k = np.arange(0, n+1)

# モデルを作成
hgeom = sp.stats.hypergeom(M, n, N)
# 確率分布を計算
pmf_hgeom = hgeom.pmf(k)

# 可視化
plt.bar(k, pmf_hgeom)
plt.xticks(k)
plt.xlabel('当たりの個数')
plt.ylabel('確率')
plt.show()
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/502109/a6421082-71f7-92ea-3cab-e291cd382a4b.png)

|当たり個数|確率|
|:---:|:---:|
|**0**|0.00010|
|**1**|0.00433|
|**2**|0.04768|
|**3**|0.19866|
|**4**|0.35759|
|**5**|0.28607|
|**6**|0.09536|
|**7**|0.01022|

##⑺ 負の二項分布（negative binomial distribution）
*   **負の二項分布**では、n回の試行でk回成功する確率を求め、最後の試行が成功であることを条件とします。
*   二項分布では成功回数kが確率変数ですが、**負の二項分布**では成功回数kは固定です。**確率変数には、試行回数nを使う代わりに失敗の数n-k**を使います。
*   例えば、コイントスをくり返して3回表が出るまでに何回投げる必要があるか。すなわち3回成功するまでに何回失敗する必要があるか、という失敗回数が確率変数になります。

```python
# パラメータを指定
N = 12  #試行回数
p = 0.5 #成功確率
k = 3   #成功回数

# 確率分布を計算
pmf_nbinom = sp.stats.nbinom.pmf(range(N), k, p)

# 可視化
plt.bar(range(N), pmf_nbinom)
plt.xlabel('失敗の回数')
plt.ylabel('確率')
plt.xticks(range(N))
plt.show()
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/502109/e30ef192-7650-4556-f4b0-d28bfaaed5ef.png)

|失敗の回数|確率|
|:---:|:---:|
|**0**|0.125|
|**1**|0.188|
|**2**|0.188|
|**3**|0.156|
|**4**|0.117|
|**5**|0.082|
|**6**|0.055|
|**7**|0.035|
|**8**|0.022|
|**9**|0.013|
|**10**|0.008|
|**11**|0.005|

##まとめ
離散型の確率分布を見てきましたが、何が確率変数となるのか、平たく言えばx軸に何を置くかという点を意識して一覧表にまとめます。

||確率分布の種類|確率変数|パラメータ|
|---|---|---|---|
|**⑴**|ベルヌーイ分布|事象0, 1|発生確率p|
|**⑵**|二項分布|試行の回数|発生確率p, 発生回数k, 試行回数n|
|**⑶**|ポアソン分布|試行の回数|平均発生回数mu|
|**⑷**|幾何分布|試行の回数|成功確率p, 試行回数k|
|**⑸**|離散一様分布|事象の種類|※scipy.atatsの一様分布は連続型のみ|
|**⑹**|超幾何分布|成功の個数|全体個数M, 全体中の成功個数n, 選択個数N|
|**⑺**|負の二項分布|失敗の回数|成功確率p, 成功回数k, 試行回数N|