LoginSignup
11
17

More than 5 years have passed since last update.

@Keisuke_Nakamura(Nakamura Keisuke)

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

Last updated at Posted at 2018-01-02

今まで統計モデリングをちゃんと勉強したことがなかったので、有名な「データ解析のための統計モデリング入門」を勉強してみようと思います。ただこの本、基本的に「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()

image.png

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

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

# 標準偏差 : 1.728040060004279

参考

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

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

  • データと重ねて表示
# 重ねて表示
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

ポアソン分布とは何か？

  • ポアソン分布
\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()

image.png

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

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

image.png

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

image.png

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

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

参考

章全体の参考

11
17
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
11
17