1. Pythonで学ぶ統計学　2. 確率分布[scipy.stats徹底理解]

*   データから計算される確率分布のことを**「経験分布」**といいます。これに対して、**確率分布を生成してくれる関数は「理論分布」**といいます。
*   まず、**分布の形（確率分布の種類）**を決める、それから**母数（確率分布のパラメータ）**を決めてしまえば、母集団分布の推定ができます。
*   そうした統計関数を集めたモジュールが**scipy.stats**です。その基本的な使い方は、次のように**記法が統一**されています。

![1_2_01.PNG](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/502109/fe16cbb3-ba6a-c5e4-c0dd-d8d4d822efd6.png)

#⑴ 確率分布の種類

*   確率関数は**「離散型」**と**「連続型」**の２つに大別されます。
*   **離散型**は、例えばサイコロの目のように**とびとびの値をとる変数**です。また**連続型**は、重量や温度のように**連続した値をとるもの**をいいます。
*   以下に、**scipy.stats**に実装されている確率分布から、知っておきたい15種類を列挙しました。

||確率分布|probability distribution|メソッド|データ|
|:---:|---|---|---|:---:|
|1|**超幾何分布**|Hypergeometric distribution|scipy.stats.**hypergeom**|<font color="DarkOrange">離散型</font>|
|2|**ベルヌーイ分布**|Bernoulli distribution|scipy.stats.**bernoulli**|<font color="DarkOrange">離散型</font>|
|3|**二項分布**|binomial distribution|scipy.stats.**binom**|<font color="DarkOrange">離散型</font>|
|4|**ポアソン分布**|Poisson distribution|scipy.stats.**poisson**|<font color="DarkOrange">離散型</font>|
|5|**幾何分布**|Geometric Distribution|scipy.stats.**geom**|<font color="DarkOrange">離散型</font>|
|6|**負の二項分布**|Negative Binomial Distribution|scipy.stats.**nbinom**|<font color="DarkOrange">離散型</font>|
|7|**一様分布**|uniform distribution|scipy.stats.**uniform**|<font color="DarkOrange">離散型</font>|
|8|**正規分布**|normal distribution|scipy.stats.**norm**|<font color="SteelBlue ">連続型</font>|
|9|**指数分布**|exponential distribution|scipy.stats.**expon**|<font color="SteelBlue ">連続型</font>|
|10|**ガンマ分布**|gamma distribution|scipy.stats.**gamma**|<font color="SteelBlue ">連続型</font>|
|11|**ベータ分布**|beta distribution|scipy.stats.**betabinom**|<font color="SteelBlue ">連続型</font>|
|12|**コーシー分布**|Cauchy distribution|scipy.stats.**cauchy**|<font color="SteelBlue ">連続型</font>|
|13|**対数正規分布**|lognormal distribution|scipy.stats.**lognorm**|<font color="SteelBlue ">連続型</font>|
|14|**パレート分布**|Pareto distribution|scipy.stats.**pareto**|<font color="SteelBlue ">連続型</font>|
|15|**ワイブル分布**|Weibull distribution|scipy.stats.**dweibull**|<font color="SteelBlue ">連続型</font>|

*   これらの確率分布は、ほぼ**共通のメソッド**を実装しています。
*   従って、メソッドを理解すれば、いろいろな確率分布を操作することが可能となります。

#⑵ 各種メソッドの機能

*   ここでは、もっとも使用頻度が高い**正規分布`scipy.stats.norm`**を例として、主要なメソッドを見ていきます。
*   引数となるのは、**入力データ`x`**、パラメータである**期待値（平均値）`loc`**と**標準偏差`scale`**のほか、**要素の数`size`**、**乱数生成のシード`random_state=整数`**などです。

```python
# 数値計算のためのライブラリ
import numpy as np
from scipy import stats

# グラフ描画のためのライブラリ
import matplotlib.pyplot as plt
%matplotlib inline

# matplotlibを日本語表示に対応させるモジュール
!pip install japanize-matplotlib
import japanize_matplotlib
```

##➀ rvs (Random variates) 確率変数
*   記法：`rvs(loc=0, scale=1, size=1, random_state=None)`
*   **確率変数**は、**確率的な法則に従って変化する値**のことです。
*   指定したパラメータに基づいて、**正規分布に従う確率変数（取り得る値）**を、**ランダムに指定した個数だけ生成**します。

```python
# rvsで正規分布に従う疑似乱数を生成
norm_rvs = stats.norm.rvs(loc=50, scale=20, size=1000, random_state=0) # 期待値=50, 標準偏差=20, 個数=1000

# 可視化（ヒストグラムに表現）
plt.hist(norm_rvs, bins=10, alpha=0.3, ec='blue')
plt.xlabel("階級", fontsize=13)
plt.ylabel("度数", fontsize=13)

plt.show()
```

![1_2_02.PNG](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/502109/227f6891-fa7a-638c-62e0-0c7ffdcf28d7.png)

##➁ pdf (Probability density function) 確率密度関数
*   記法：`pdf(x, loc=0, scale=1)`
*   **確率密度**は、定義された域内での確率変数Xの値の**相対的な出やすさ**を表します。
*   平たく言えば、**確率密度関数は、<font color="OrangeRed">連続型</font>のデータを引数にとると確率密度が算出される関数**のことです。

```python
# 等差数列を生成
X = np.arange(start=1, stop=7.1, step=0.1)

# pdfで確率密度関数を生成
norm_pdf = stats.norm.pdf(x=X, loc=4, scale=0.8) # 期待値=4, 標準偏差=0.8 

# 可視化
plt.plot(X, norm_pdf)
plt.xlabel("確率変数X", fontsize=13)
plt.ylabel("確率密度pdf", fontsize=13)
plt.show()
```

![1_2_03.PNG](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/502109/eaa12393-385e-79d8-db9e-bc7fb33e1821.png)

*   試みに、確率変数 x=5 のときの確率密度を確認してみます。

```python
# 確率変数5のときの確率密度を計算
x = 5
y = stats.norm.pdf(x=x, loc=4, scale=0.8) # 期待値=4, 標準偏差=0.8
print("確率変数x=5のときの確率密度：", y)

# 確率密度関数のグラフにプロット
plt.plot(X, norm_pdf)
plt.plot(x, y, 'bo') # 青色ドットを布置

plt.vlines(x, 0.0, y, lw=2, linestyles='dashed') # 垂直線
plt.hlines(y, 1.0, x, lw=2, linestyles='dashed') # 水平線

plt.xlabel("確率変数X", fontsize=13)
plt.ylabel("確率密度pdf", fontsize=13)

plt.show()
```

![1_2_04.PNG](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/502109/b397202e-ae45-742f-9af0-9939997eef86.png)

*   **<font color="OrangeRed">確率密度関数pdfは連続型データのためのメソッドであり、離散型データには次の確率質量関数pmfを使います。</font>**

##➂ pmf (Probability mass function) 確率質量関数
*   記法：`pmf(k, n, p, loc=0)`
*   **確率質量**は、確率変数Xのとびとびの要素ごとの**相対的な出やすさ**を表します。
*   平たく言えば、**確率質量関数は、<font color="OrangeRed">離散型</font>のデータを引数にとると確率質量が算出される関数**のことです。
*   ここでは**２項分布**を例として、「成功確率40％の試行を５回行ったとき、そのうち成功する回数（０～４回）ごとの確率」を示します。

```python
# 成功確率
p = 0.4
# 試行回数
n = 5
# 成功回数
k = np.arange(0, 5) # array([0, 1, 2, 3, 4])

# pmfで成功回数ごとの確率を計算
binom_pmf = stats.binom.pmf(k, n, p) # 成功回数=0～4, 試行回数=5, 成功確率=0.4

# 可視化
plt.plot(k, binom_pmf, 'bo', ms=8)
plt.vlines(k, 0, binom_pmf, colors='b', lw=3, alpha=0.5)

plt.xticks(k) # x軸目盛

plt.xlabel("確率変数X", fontsize=13)
plt.ylabel("確率質量関数pmf", fontsize=13)

plt.show()
```

![1_2_14.PNG](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/502109/5a98870a-16df-29e5-f17d-e008728e6c2f.png)

##➃ logpdf (Log of the probability density function) ログ確率密度関数
*   記法：`logpdf(x, loc=0, scale=1)`
*   **ログ確率密度関数**は、**底をeとする確率密度の対数**をとったものです。
*   ネイピア数 $e=2.7182818284 ...$ と続く超越数を$a$乗したら確率密度になるときの$a$の値ことです。式に表すと $e^a=$pdf です。

```python
# 等差数列を生成
X = np.arange(start=1, stop=7.1, step=0.1)

# logpdfで確率密度関数の対数を生成
norm_logpdf = stats.norm.logpdf(x=X, loc=4, scale=0.8) # 期待値=4, 標準偏差=0.8

# 可視化
plt.plot(X, norm_logpdf)
plt.xlabel("確率変数X", fontsize=13)
plt.ylabel("ログ確率密度関数logpdf", fontsize=13)
plt.show()
```

![1_2_05.PNG](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/502109/7fed1c5c-236f-f2e3-79b1-abab85fbb058.png)

*   試みに、上図で確率密度の対数が最大値を示す確率変数 x=4 について確認してみます。
*   手順としては、まず確率変数 **x=4 のログ確率密度の値**を取得し、次いで **x=4 の確率密度を計算してからその対数**をとり、両方を比較します。

```python
# ログ確率密度の最大値を取得
logpdf_max = np.max(norm_logpdf) 
print('x=4のログ確率密度：', logpdf_max)

# 確率密度を計算
x = 4
v_pdf = stats.norm.pdf(x=x, loc=4, scale=0.8) # 期待値=4, 標準偏差=0.8

# 確率密度の対数を計算
v_pdf_log = np.log(v_pdf)
print('x=4の確率密度の対数：', v_pdf_log)

# 可視化
plt.plot(X, norm_logpdf)
plt.plot(4, logpdf_max, 'bo') # 青色ドットを布置

logpdf_min = np.min(norm_logpdf) # ログ確率密度の最小値(描画用)
plt.vlines(x, logpdf_max, logpdf_min, lw=2, linestyles='dashed') # 垂直線
plt.hlines(logpdf_max, 1.0, x, lw=2, linestyles='dashed') # 水平線

plt.xlabel("確率変数X", fontsize=13)
plt.ylabel("ログ確率密度関数logpdf", fontsize=13)

plt.show()
```

![1_2_06.PNG](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/502109/0ae54a1c-0aa0-6129-233f-d07ea43cb941.png)

##➄ cdf (Cumulative distribution function) 累積分布関数
*   記法：`cdf(x, loc=0, scale=1)`
*   **累積分布関数**は、**確率変数$X$がある値$x$以下となる確率**を計算します。
*   例えば、サイコロを投げたときに「出る目が３以下」となる確率は、１から３までの確率密度をすべて足し合わせたものになります。

```python
# x軸の等差数列を生成
X = np.linspace(start=1, stop=7, num=100)

# cdfで累積分布関数を生成
norm_cdf = stats.norm.cdf(x=X, loc=4, scale=0.8) # 期待値=4, 標準偏差=0.8

# cdfでx=5以下の累積確率を計算
under_5 = stats.norm.cdf(x=5, loc=4, scale=0.8) # 期待値=4, 標準偏差=0.8
print('5以下になる確率：', under_5)

# 可視化
plt.plot(X, norm_cdf)
plt.plot(5, under_5, 'bo') # 青色ドットを布置

plt.vlines(5, 0.0, under_5, lw=2, linestyles='dashed') # 垂直線
plt.hlines(under_5, 1.0, 5, lw=2, linestyles='dashed') # 水平線

plt.xlabel("確率変数X", fontsize=13)
plt.ylabel("累積確率密度cdf", fontsize=13)

plt.show()
```

![1_2_07.PNG](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/502109/bb0a3c8a-d282-1968-012e-a80e7d97ff16.png)

*   全く同一のパラメータによる確率密度関数の上に、同じく確率変数 x=5 のポイントを示したのが下図になります。
*   累積分布関数は、この図の水色の領域に相当します。

![1_2_08.PNG](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/502109/6e5009d7-4c51-0036-2652-a085e0124306.png)

##➅ sf (Survival function) 生存関数

*   記法：`sf(x, loc=0, scale=1)`
*   **生存関数**は、累積分布関数cdfとは逆に、**確率変数$X$がある値$x$<font color="OrangeRed">以上</font>となる確率**を計算します。
*   **｢1-cdf｣** と定義されることもありますが、**sfの方がより正確な場合もある**との指摘がscipy.org https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.norm.html に見られます。

![1_2_15.PNG](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/502109/cd8b3e62-b215-5178-a75c-9614ac738447.png)

```python
# x軸の等差数列を生成
X = np.linspace(start=1, stop=7, num=100)

# sfで生存関数を生成
norm_sf = stats.norm.sf(x=X, loc=4, scale=0.8) # 期待値=4, 標準偏差=0.8

# sfでx=5以上の累積確率を計算
upper_5 = stats.norm.sf(x=5, loc=4, scale=0.8) # 期待値=4, 標準偏差=0.8
print('5以上になる確率：', upper_5)

# 可視化
plt.plot(X, norm_sf)
plt.plot(5, upper_5, 'bo') # 青色ドットを布置

plt.vlines(5, 0.0, upper_5, lw=2, linestyles='dashed') # 垂直線
plt.hlines(upper_5, 1.0, 5, lw=2, linestyles='dashed') # 水平線

plt.xlabel("確率変数X", fontsize=13)
plt.ylabel("生存関数sf", fontsize=13)

plt.show()
```

![1_2_09.PNG](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/502109/773e3813-c3e5-9359-3f53-46f52bd8bc24.png)

*   全く同一のパラメータによる確率密度関数の上に、同じく確率変数 x=5 のポイントを示したのが下図です。
*   生存関数は、図中の水色の領域に相当します。

![1_2_10.PNG](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/502109/cd5a29de-ab4e-cf81-aae2-a6a94546dc4e.png)

##➆ ppf (Percent point function) パーセント点関数

*   記法：`ppf(q, loc=0, scale=1)`
*   **パーセント点関数**は、累積分布関数cdfの逆関数で、**累積分布関数をq％と指定するとその値をとる変数**を返します。
*   従って、ppf(0.25)は第１四分位点、ppf(0.5)は中央値にあたる第２四分位点、ppf(0.75)は第３四分位点に相当します。

```python
# x軸の等差数列を生成
X = np.arange(start=1, stop=7, step=0.1)

# pdfで確率密度関数を生成
norm_pdf = stats.norm.pdf(x=X, loc=4, scale=0.8) # 期待値=4, 標準偏差=0.8

# ppfで累積分布関数75％に当たる変数を取得
q_75 = stats.norm.ppf(q=0.75, loc=4, scale=0.8) # 期待値=4, 標準偏差=0.8
print('累積分布関数75％点の確率変数：', q_75)

# pdfで当該変数の確率密度を取得
v = stats.norm.pdf(x=q_75, loc=4, scale=0.8)
print('累積分布関数75％点の確率密度：', v)

# 可視化
plt.plot(X, norm_pdf)
plt.plot(q_75, v, 'bo') # 青色ドットを布置

plt.vlines(q_75, 0.0, v, lw=2, linestyles='dashed') # 垂直線
plt.hlines(v, 1.0, q_75, lw=2, linestyles='dashed') # 水平線

plt.xlabel("確率変数X", fontsize=13)
plt.ylabel("確率密度関数pdf", fontsize=13)

plt.show()
```

![1_2_11.PNG](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/502109/238e8805-bd7b-ee36-4a0f-bc7946071923.png)

##➇ isf (Inverse survival function) 逆生存関数

*   記法：`isf(q, loc=0, scale=1)`
*   **逆生存関数**は、生存関数sfの逆関数ですが、やることはパーセント点関数ppfと逆のことです。**生存関数の残余をq％と指定するとその値をとる変数**を返します。
*   従って、四分位点はパーセント点関数ppfとは逆になり、isf(0.25)は第３四分位点、isf(0.5)は同じく第２四分位点、isf(0.75)は第１四分位点に相当します。

```python
# x軸の等差数列を生成
X = np.arange(start=1, stop=7, step=0.1)

# pdfで確率密度関数を生成
norm_pdf = stats.norm.pdf(x=X, loc=4, scale=0.8) # 期待値=4, 標準偏差=0.8

# isfで逆生存関数25％に当たる変数を取得
q_25 = stats.norm.isf(q=0.25, loc=4, scale=0.8) # 期待値=4, 標準偏差=0.8
print('逆生存関数25％点の確率変数：', q_25)

# pdfで当該変数の確率密度を取得
v = stats.norm.pdf(x=q_25, loc=4, scale=0.8)
print('逆生存関数25％点の確率密度：', v)

# 可視化
plt.plot(X, norm_pdf)
plt.plot(q_25, v, 'bo') # 青色ドットを布置

plt.vlines(q_25, 0.0, v, lw=2, linestyles='dashed') # 垂直線
plt.hlines(v, 1.0, q_25, lw=2, linestyles='dashed') # 水平線

plt.xlabel("確率変数X", fontsize=13)
plt.ylabel("確率密度関数pdf", fontsize=13)

plt.show()
```

![1_2_12.PNG](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/502109/a7774926-4047-9de9-bfb9-5b06db3441c7.png)

##➈ interval (Endpoints of the range that contains alpha percent of the distribution) 区間推定

*   記法：`interval(alpha, loc=0, scale=1)`
*   実際のデータから計算される**経験分布**では、平均値も標準偏差も常に未知のものです。標本から得られた平均値が、母集団における**真の母数$θ$**に近似できているか本当はわかりません。
*   そこで、**一定の確率のもとに<font color="OrangeRed">母数$θ$を含む区間</font>を求めることを区間推定**といいます。

```python
# x軸の等差数列を生成
X = np.arange(start=1, stop=7, step=0.1)

# pdfで確率密度関数を生成
norm_pdf = stats.norm.pdf(x=X, loc=4, scale=0.8) # 期待値=4, 標準偏差=0.8

# intervalで信頼区間95％に当たる変数を取得
lower, upper = stats.norm.interval(alpha=0.95, loc=4, scale=0.8) # 信頼係数=0.95, 期待値=4, 標準偏差=0.8
print('信頼区間95％の下限：', lower)
print('信頼区間95％の上限：', upper)

# pdfで各変数の確率密度を取得
v_lower = stats.norm.pdf(x=lower, loc=4, scale=0.8) # 期待値=4, 標準偏差=0.8
v_upper = stats.norm.pdf(x=upper, loc=4, scale=0.8) # 期待値=4, 標準偏差=0.8

# 可視化
plt.plot(X, norm_pdf)
plt.plot(lower, 0.0, 'k|')
plt.plot(upper, 0.0, 'k|')

plt.vlines(lower, 0.0, v_lower, lw=0.8) # 下限の垂直線
plt.vlines(upper, 0.0, v_upper, lw=0.8) # 上限の垂直線

plt.xlabel("確率変数X", fontsize=13)
plt.ylabel("確率密度関数pdf", fontsize=13)

plt.show()
```

![1_2_13.PNG](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/502109/b5de5e62-98a1-2b1a-a1d1-bedb1092b7ea.png)

*   区間推定では、**まず母数$θ$を含む区間の確率を指定**します。この例では`alpha=0.95`がそれで、これを**信頼係数**と呼びます。多くの場合、99％、95％、90％が使われますが、これは任意です。
*   その指定された確率のもとで**母数$θ$を含む区間の下限と上限**が求められます。この例では95％の確率で、母数$θ$は確率変数が2.43～5.57の間に含まれていることを示しています。


#あとがき
*   ログ累積分布関数logcdf、ログ生存関数logsfなどは省きましたが、およそ使用頻度の高そうなものは押さえられたと思います。
*   これらを利用して、先に掲げた15種類の確率分布を一つずつ解きほぐしてみたいと思っています。
*   [1. Pythonで学ぶ統計学　2-1. 確率分布［離散型変数］](https://qiita.com/y_itoh/items/cf1229bebc4839041fd5)