Pythonで確率プロット

# Pythonで確率プロット

## 目的

よく業務で使用しているPythonを用いた確率プロットの作成方法について紹介したい。
(検索しても良い情報がなかったので…)

参考: [正規確率プロット（Normal Q-Q Plot）の作成 with Excel](https://hitorimarketing.net/tools/normal-quantile-quantile-plot.html)

こちらはExcelで確率プロットを作成しているページ。
本記事で紹介しているのは、ほぼこちらのPython版となる。

## 確率プロットとは

QQ-Plotとも言う。
簡単に言えば、データが正規分布に即しているかどうかを視覚的に判断できるようにしたもの。
正規分布から外れているデータを発見したり、近似の直線を引くことで、データがどのくらいの範囲に収まるかを予測したりする。

正規分布に関しては以下の書籍がわかりやすくて勉強になる。

<div class="amazlet-box" style="margin-bottom:0px;"><div class="amazlet-image" style="float:left;margin:0px 12px 1px 0px;"><a href="http://www.amazon.co.jp/exec/obidos/ASIN/4062579855" name="amazletlink" target="_blank"><img src="https://images-fe.ssl-images-amazon.com/images/I/51W7Xg0yoiL._SL160_.jpg" alt="経済数学の直観的方法 確率・統計編 (ブルーバックス)" style="border: none;" /></a></div><div class="amazlet-info" style="line-height:120%; margin-bottom: 10px"><div class="amazlet-name" style="margin-bottom:10px;line-height:120%"><a href="http://www.amazon.co.jp/exec/obidos/ASIN/4062579855" name="amazletlink" target="_blank">経済数学の直観的方法 確率・統計編 (ブルーバックス)</a><div class="amazlet-powered-date" style="font-size:80%;margin-top:5px;line-height:120%">posted with <a href="http://www.amazlet.com/" title="amazlet" target="_blank">amazlet</a> at 18.12.16</div></div><div class="amazlet-detail">長沼 伸一郎 <br />講談社 <br />売り上げランキング: 82,727<br /></div><div class="amazlet-sub-info" style="float: left;"><div class="amazlet-link" style="margin-top: 5px"><a href="http://www.amazon.co.jp/exec/obidos/ASIN/4062579855" name="amazletlink" target="_blank">Amazon.co.jpで詳細を見る</a></div></div></div><div class="amazlet-footer" style="clear: left"></div></div>

## 環境について

GoogleColaboratryのPython3ノートブックで作成した。
JupyterNotebookでも使えると思う。

## コード


### 使用するライブラリをimport

```python
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy import stats
```

### ndarrayとDataFrameでデータを用意する

今回はサンプル用に正規分布に則ったデータを作成する。実際の業務ではcsvファイルから読み込んだりすることになる。
使用しやすいようにndarrayとdataframeの両方で用意する。
また、列数と行数を確保しておく。

```python
# 正規分布のデータを3つ作成
npspl = np.random.normal(50, 10, [500, 3])
rownum, colnum = npspl.shape[0], npspl.shape[1]
dfspl = pd.DataFrame(npspl, columns=['Data1','Data2','Data3'])
```
### ヒストグラムプロット

データが正規分布であることを確認するためヒストグラムをプロットしてみる。

```python
for i in range(colnum):
  plt.hist(npspl[:,i], 20)
  plt.show()
```

以下のように出力される。

![ダウンロード.png](https://qiita-image-store.s3.amazonaws.com/0/277413/3178a8cc-86c2-c58b-1651-8ad3eacbe0f0.png)

### パーセンタイル順位を作成、Normsinvに変換する

ここが確率プロットのポイント。

* **パーセンタイル順位とは？**

順位をパーセントに変換したもの。
例えば、とある値が全体の上位20%の大きさだったらパーセンタイル順位は0.8となる。

* **Normsinvとは？**

標準正規分布における期待値で、パーセンタイル順位から算出する。


```python
norspl = stats.norm.ppf((
    np.array(dfspl.rank(method='first', na_option='keep') - 1/2)/ rownum),
    loc=0, scale=1)
```

ワンライナーで書いちゃったので分解して解説する。

以下が各値のパーセンタイル順位を算出している部分。

```python
np.array(dfspl.rank(method='first', na_option='keep') - 1/2)/ rownum)
```

最初にDataFrameのrankを用いて各値の順位を求める。(参考: [DatarFrameのrankのリファレンス](https://pandas.pydata.org/pandas-docs/stable/generated/pandas.DataFrame.rank.html))

その後、扱いやすいようndarrayに変換した後、データ数で割ってパーセンタイル順位をはじき出している。

-1/2しているのは、R言語のqqplot関数を用いるとこのような計算をしているから。
個人的には少し謎な部分で、-1/2しようがしまいが特に変わらないと思う。

上の式の外側の部分が以下。

```python
stats.norm.ppf('上で記載した式', loc=0, scale=1)
```

ここでパーセンタイル順位がそれぞれ標準正規分布のどこに位置するかどうかの期待値を算出している。(参考: [norm.ppfのリファレンス](https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.norm.html))

loc=0が平均の設定で、scale=1が分散の設定。

### グラフをプロット

#### 各種設定

色やグラフのサイズを設定する。

```python
colorsets=['b','g','r']
xmin, xmax = -6, 6
ymin, ymax = 0, 100
xgrid, ygrid = 2, 10
```

#### 散布図でプロット

```python
for i in range(colnum):
  dataname = dfspl.columns[i]
  plt.scatter(norspl[:,i], npspl[:,i], color = colorsets[i], label = dataname)
plt.xlim([xmin, xmax])
plt.xticks(np.arange(xmin, xmax+1, xgrid))
plt.ylim([ymin, ymax])
plt.yticks(np.arange(ymin, ymax+1, ygrid))
plt.title('QQ-Plot')
plt.grid(linestyle='dashed')
plt.legend(loc=2)
plt.show()
```

### 出力された確率プロット

以下のように線形になっていれば、正規分布に即したデータだと判断できる。

![ダウンロード (1).png](https://qiita-image-store.s3.amazonaws.com/0/277413/b76890a2-0ba2-739f-afe1-a91b11fca7e3.png)