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

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

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

データの取得

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

# サンプルデータの取得
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として表示)

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]

データ数を調べる

print(len(data))
#50

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

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

度数分布を取得

data.value_counts(sort=False)

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

ヒストグラムを作成

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()

標本分散

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

#分散 : 2.986122448979592

標本標準偏差

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

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

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()

データと重ねて表示

# 重ねて表示
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()

ポアソン分布とは何か？

ポアソン分布

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

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

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()

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

図2.7

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()

図2.8

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()

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

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()

参考

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