データ解析のための統計モデリング入門(for Python) : ２章

今まで統計モデリングをちゃんと勉強したことがなかったので、有名な「データ解析のための統計モデリング入門」を勉強してみようと思います。ただこの本、基本的に「R」を使うことを前提に書かれてるわけですが、個人的にPythonが好きなので、出来る限りPythonを使って再現していこうと思います。もちろん同じことをすでにやっている人たちはたくさんいるわけですが、自分でまとめることに意味があると思うので、とりあえずやっていこうと思います。

# 例題 : 粒子数の統計モデリング

サンプルデータの内容を把握します。必要なデータはhttp://hosho.ees.hokudai.ac.jp/~kubo/ce/IwanamiBook.html よりダウンロードし適当な場所においておきます。

- データの取得

Rのデータのため、pypeRを使ってRコマンドを実行し、それをpython側にpandas形式で取り込みます

```python3
# サンプルデータの取得
import pandas as pd
import pyper as pr

workPath = "サンプルデータをおいた場所"
r = pr.R(use_numpy='True', use_pandas='True') #インスタンス作成
r('load("{}/data.RData")'.format(workPath)) #実行
data = pd.Series(r.get('data')) #pandasのSeriesとして格納
```

- データの把握をする(numpyのarrayとして表示)

```python3
import numpy as np
dataNP = np.array(data)
print(dataNP)

#[2 2 4 6 4 5 2 3 1 2 0 4 3 3 3 3 4 2 7 2 4 3 3 3 4 3 7 5 3 1 7 6 4 6 5 2 4
# 7 2 2 6 2 4 5 4 5 1 3 2 3]
```

- データ数を調べる

```python3
print(len(data))
#50
```

- データの要約(jupyter notebookで表示)

```python3
data.describe()

#count    50.00000
#mean      3.56000
#std       1.72804
#min       0.00000
#25%       2.00000
#50%       3.00000
#75%       4.75000
#max       7.00000
#dtype: float64
```

- 度数分布を取得

```python3
data.value_counts(sort=False)

#0     1
#1     3
#2    11
#3    12
#4    10
#5     5
#6     4
#7     4
#dtype: int64
```

- ヒストグラムを作成

```python3
import matplotlib.pyplot as plt
plt.hist(data, bins=np.arange(-0.5, 8.5, 1.0), ec="black")
plt.xlabel('data', fontsize=12)
plt.ylabel('Frequency', fontsize=12)
plt.title('Histgram of data', fontsize=12)
plt.xticks([2*x for x in range(5)])
plt.show()
```
![image.png](https://qiita-image-store.s3.amazonaws.com/0/216851/80a866a1-107f-9ca0-f002-adebb05b0917.png)

- 標本分散

```python3
print("分散 : {}".format(data.var()))

#分散 : 2.986122448979592
```

- 標本標準偏差

```python3
print("標準偏差 : {}".format(np.sqrt(data.var())))
# print("標準偏差 : {}".format(np.std(data, ddof=1))) こちらでもOK

# 標準偏差 : 1.728040060004279
```

#### 参考
- pythonのnumpy.arrayの統計量一覧を表示するワンライナー : https://qiita.com/AnchorBlues/items/051dc69e81705b52adad
- ヒストグラム : https://stats.biopapyrus.jp/python/hist.html

# データと確率分布の対応関係をながめる

```python3
from scipy import stats
import matplotlib.font_manager
# 日本語表示を可能にする
fontprop = matplotlib.font_manager.FontProperties(fname="/Library/Fonts/ヒラギノ丸ゴ ProN W4.ttc")
# 粒子数がyである確率をリストとして保存
y = np.arange(0, 10)
prob = pd.Series(stats.poisson.pmf(y, 3.56), index=y)
# 折れ線として表示
plt.plot(prob, 'o--')
plt.xlabel('$y$', fontsize=12)
plt.ylabel('prop', fontsize=12)
plt.title('平均$\lambda=3.56$のポアソン分布', fontdict={'fontproperties':fontprop}, fontsize=12)
plt.show()
```
![image.png](https://qiita-image-store.s3.amazonaws.com/0/216851/5f85ad0f-12c7-9e9b-79f1-b9dcbc64f0cd.png)

- データと重ねて表示

```python3
# 重ねて表示
plt.plot(prob * 50, 'o--')
plt.hist(data, bins=np.arange(-0.5, 8.5, 1.0), ec="black")
plt.xlabel('data', fontsize=12)
plt.ylabel('Frequency', fontsize=12)
plt.title('平均$\lambda=3.56$のポアソン分布', fontdict={'fontproperties':fontprop}, fontsize=12)
plt.show()
```
![image.png](https://qiita-image-store.s3.amazonaws.com/0/216851/cb5601b3-2c71-244b-2818-6463ee0ac21d.png)

# ポアソン分布とは何か？

- ポアソン分布

```math
\begin{aligned}
p(y\ |\ \lambda)=\frac{\lambda^{y}\exp{(-\lambda)}}{y!}
\end{aligned}
```

- 平均$\lambda$を変えたときのポアソン分布

```python3
import matplotlib.ticker as tick
y1 = range(21)
prob1 = pd.Series(stats.poisson.pmf(y1, 3.5), index=y1)
prob2 = pd.Series(stats.poisson.pmf(y1, 7.7), index=y1)
prob3 = pd.Series(stats.poisson.pmf(y1, 15.1), index=y1)
matplotlib.pyplot.rcParams['figure.figsize'] = (7.0, 5.0)
fontprop = matplotlib.font_manager.FontProperties(fname="/Library/Fonts/ヒラギノ丸ゴ ProN W4.ttc")
plt.title("さまざまな平均($\lambda$)のポアソン分布", fontdict = {"fontproperties": fontprop}, fontsize=12)
plt.plot(prob1, 'o--')
plt.plot(prob2, '^--')
plt.plot(prob3, 'd--')
plt.xlabel("y", fontsize=12)
plt.ylabel("prob", fontsize=12)
plt.legend(['3.5','7.7', '15.1'])
plt.xticks([0, 5, 10, 15, 20])
plt.yticks([0.00, 0.10, 0.20])
plt.gca().yaxis.set_minor_locator(tick.MultipleLocator(0.05))
plt.show()
```

![image.png](https://qiita-image-store.s3.amazonaws.com/0/216851/6897785b-8c07-6071-64b7-56248c6cbc99.png)

# ポアソン分布のパラメータの最尤推定
- 図2.7

```python3
matplotlib.pyplot.rcParams['figure.figsize'] = (15.0, 12.0)
x = np.arange(2.0, 5.6, 0.4)
y = np.arange(0, 10)
for i, _x in enumerate(x):
    plt.subplot(3,3,i+1)
    prob = pd.Series(stats.poisson.pmf(y, _x), index=y)
    plt.plot(prob * 50, 'o--')
    plt.hist(data, bins=np.arange(-0.5, 8.5, 1.0), ec="black")
    plt.yticks([0,5,10,15])
    plt.legend(['lambda={}'.format(round(_x,1)), 'logL={}'.format(round(sum(stats.poisson.logpmf(data, _x))), 1)])
plt.show()
```
![image.png](https://qiita-image-store.s3.amazonaws.com/0/216851/223ae17b-f3a4-f347-2580-5e33b6e6352b.png)

- 図2.8

```python3
fontprop = matplotlib.font_manager.FontProperties(fname="/Library/Fonts/ヒラギノ丸ゴ ProN W4.ttc")
matplotlib.pyplot.rcParams['figure.figsize'] = (7.0, 5.0)
x = np.arange(2.0, 5.1, 0.1)
y = [sum(stats.poisson.logpmf(data, _x)) for _x in x]
varlam = pd.Series(y, index=x)

plt.plot(varlam, '-')
plt.xticks([n for n in np.arange(2.0, 5.5, 0.5)])
plt.yticks([-120, -110, -100])
plt.xlim([1.9,5.1])
plt.ylim([-123, -96])
plt.vlines([3.56],-125, -95,  "blue", linestyles='dashed') 
plt.xlabel("$\lambda$", fontsize=12,)
plt.ylabel("対数尤度", fontdict = {"fontproperties": fontprop}, fontsize=12)
plt.gca().yaxis.set_minor_locator(tick.MultipleLocator(5.0))
plt.show()
```
![image.png](https://qiita-image-store.s3.amazonaws.com/0/216851/c28a546e-2f36-ac58-8ede-3ce70efa6c35.png)

# 擬似乱数と最尤推定値のばらつき

```python3
fontprop = matplotlib.font_manager.FontProperties(fname="/Library/Fonts/ヒラギノ丸ゴ ProN W4.ttc")
matplotlib.pyplot.rcParams['figure.figsize'] = (7.0, 5.0)

result = []
for i in range(3000):
    dummyData = np.random.poisson(lam=3.5, size=50)
    x_sum = np.sum(dummyData)
    varlam = 1/50 * x_sum
    result.append(varlam)

result = pd.Series(result)
plt.hist(result, bins=20, ec='black')
plt.xticks([n for n in np.arange(2.5, 5.0, 0.5)])
plt.yticks([0,100,300])
plt.gca().yaxis.set_minor_locator(tick.MultipleLocator(100))
plt.xlabel('試行ごとに推定された$\lambda$', fontdict={"fontproperties":fontprop}, fontsize=12)
plt.title(r'ポアソン分布の平均$\lambda$の最尤推定値$\bar{\lambda}$のばらつき', fontdict={'fontproperties':fontprop}, fontsize=12)
plt.show()
```

![image.png](https://qiita-image-store.s3.amazonaws.com/0/216851/f70b74c9-d240-e1c1-315a-38ed984b996a.png)

#### 参考
- Numpyによる乱数生成まとめ : https://qiita.com/yubais/items/bf9ce0a8fefdcc0b0c97

# 章全体の参考
- 生態学のデータ解析 - 本/データ解析のための統計モデリング入門 : http://hosho.ees.hokudai.ac.jp/~kubo/ce/IwanamiBook.html
- Pythonで「データ解析のための統計モデリング入門」：http://imaimamu.com/archives/1875
- データ解析のための統計モデリング入門 : https://qiita.com/makora9143/items/6ec424776855cfafd724
- matplotlib入門 グラフの体裁を整える : http://yubais.net/doc/matplotlib/modify.html
- matplotlib入門 とりあえず描く : http://yubais.net/doc/matplotlib/introduction.html
- タオルケット体操 覚えるだけでPythonのコードが少し綺麗になる頻出イディオム : http://hachibeechan.hateblo.jp/entry/Python-idiom-101
- ipython notebookでmatplotlibでプロットした図が小さい : https://github.com/eisoku9618/kuroiwa_demos/issues/69