「データ解析のための統計モデリング入門」、もっと早くに取り組んでおけば良かった。データ分析業のアプローチに必要な観点を学べる良テキストだ。しかし本文中はRが使われているのでPythonで一通り書き直して読み進めてみる。
GLMフィッティングの際の勾配降下、確率的勾配降下法とかいろいろあったと記憶しているが、最もナイーブな実装しか書けなかった。修行が足りない。しかしPythonのstatsmodelsは今回初めて使ったが便利すぎる。
GLMフィッティングの際の勾配降下、確率的勾配降下法とかいろいろあったと記憶しているが、最もナイーブな実装しか書けなかった。修行が足りない。しかしPythonのstatsmodelsは今回初めて使ったが便利すぎる。
 | データ解析のための統計モデリング入門――一般化線形モデル・階層ベイズモデル・MCMC (確率と情報の科学) 久保 拓弥 岩波書店 売り上げランキング : 3193 Amazonで詳しく見る by AZlink |
3.3 統計モデリングの前にデータを図示する¶
In [5]:
df = pd.read_csv('./data3a.csv')
In [12]:
def draw_scatter(df):
plt.scatter(df[df.f=='C'].x, df[df.f=='C'].y, color='r')
plt.scatter(df[df.f=='T'].x, df[df.f=='T'].y, color='b')
plt.legend([u'f == C', u'f == T(施肥)'], loc='upper left')
In [13]:
draw_scatter(df)
plt.title(u'図3.2')
Out[13]:
In [15]:
df.groupby(df.f).boxplot(column='y')
plt.ylim(0, 15)
Out[15]:
3.4.1 線形予測と対数リンク関数¶
In [48]:
X = np.linspace(-4, 5, 100)
plt.plot(X, np.exp(-2 - 0.8*X), 'k--')
plt.plot(X, np.exp(-1 + 0.4*X), 'k-')
plt.legend(['$\{beta_1, beta_2\} = \{-2, -0.8\}$', '$\{beta_1, beta_2\} = \{-1, 0.4\}$'])
plt.title(u'図3.4')
Out[48]:
3.4.2 あてはめと、あてはまりの良さ¶
In [49]:
def log_likelifood(b1, b2, x, y):
# ポワソン回帰の対数尤度を返す
ret = 0
for i in range(x.size):
ret += y[i] * (b1 + b2*x[i]) - np.exp(b1 + b2*x[i]) - sum(np.log(range(1, y[i])))
return ret
独自で対数尤度最大化を実装してみる
In [59]:
def climb_b1(b1, b2, x, y):
# b2固定でb1を動かした時の最大の対数尤度を返す
sign = +1
alpha = 0.01
logL = log_likelifood(b1, b2, x, y)
while True:
b1 += alpha * sign
logL_n = log_likelifood(b1, b2, x, y)
if np.abs(logL - logL_n) < 0.00001:
break
if logL_n < logL:
sign = sign * -1
alpha = alpha * 0.5
else:
alpha = alpha * 1.2
logL = logL_n
return logL, b1
def hill_climb(alpha, x, y):
b1 = 1
b2 = 0.1
sign = +1
max_logL, b1 = climb_b1(b1, b2, x, y)
while True:
b2 += alpha * sign
# b2を動かした時の対数尤度を出す
max_logL_n, b1 = climb_b1(b1, b2, x, y)
# 変化しなくなったら終了
if np.abs(max_logL_n - max_logL) < 0.00001:
break
if max_logL_n < max_logL:
# 対数尤度が減ったら符号逆転
sign = sign * -1
alpha = alpha * 0.5
else:
alpha = alpha * 1.2
max_logL = max_logL_n
print('Likelihood: %s, b1: %s, b2: %s' % (max_logL, b1, b2))
return b1, b2
hill_climb(0.1, df.x.values, df.y.values)
Out[59]:
statsmodelsを使うパターン
In [60]:
import statsmodels.api as sm
from statsmodels.formula.api import glm
def calc_x_related_model(df):
return glm('y ~ x', data=df, family=sm.families.Poisson(sm.families.links.log)).fit()
In [61]:
x_related_model = calc_x_related_model(df)
x_related_model.summary()
Out[61]:
In [62]:
import scipy.stats
def plot_fig36(model):
X = np.linspace(-0.1, 1.6, 200)
plt.plot(X, scipy.stats.norm.pdf(X, loc=model.params[1], scale=model.bse[1]))
plt.plot(X, scipy.stats.norm.pdf(X, loc=model.params.Intercept, scale=model.bse.Intercept))
plt.title(u'図3.6')
plot_fig36(x_related_model)
In [63]:
def plot_likelifood():
X = np.linspace(0.05, 0.10, 1000)
Y = []
for x in X:
Y.append(log_likelifood(1.2926, x, df.x.values, df.y.values))
plt.plot(X, Y)
plot_likelifood()
plt.title(u'尤度が最大になる付近のグラフは正規分布に近い?')
plt.xlabel('b1')
plt.ylabel(u'対数尤度')
Out[63]:
対数尤度の評価
In [64]:
x_related_model.llf
Out[64]:
AICで評価
In [65]:
x_related_model.aic
Out[65]:
3.4.3 ポアソン回帰モデルによる予測¶
In [66]:
def draw_estimate_x(model):
X = np.linspace(7, 13, 100)
Y = model.predict({'x': X})
plt.plot(X, Y, 'k-')
In [67]:
draw_scatter(df)
draw_estimate_x(x_related_model)
plt.xlim(7, 13)
plt.title(u'体サイズ$xi$を使った統計モデルによる予測')
Out[67]:
3.5 説明変数が因子型の統計モデル¶
In [68]:
def calc_f_related_model(df):
return glm('y ~ f', data=df, family=sm.families.Poisson(sm.families.links.log)).fit()
In [69]:
f_related_model = calc_f_related_model(df)
f_related_model.summary()
Out[69]:
In [70]:
def draw_estimate_f(model):
X = np.linspace(7, 14, 1000)
Y = model.predict({'x': X, 'f': ['C']*len(X)})
plt.plot(X, Y, 'r-')
Y = model.predict({'x': X, 'f': ['T']*len(X)})
plt.plot(X, Y, 'b-')
In [71]:
draw_scatter(df)
draw_estimate_f(f_related_model)
plt.xlim(7, 14)
plt.ylim(4, 10)
plt.title(u'施肥効果$fi$を使った統計モデルによる推定')
Out[71]:
3.6 説明変数が数量型 + 因子型の統計モデル¶
In [72]:
def calc_x_f_related_model(df):
return glm('y ~ x + f', data=df, family=sm.families.Poisson(sm.families.links.log)).fit()
In [73]:
x_f_related_model = calc_x_f_related_model(df)
x_f_related_model.summary()
Out[73]:
In [74]:
def draw_estimate_x_f(model):
X = np.linspace(7, 13, 100)
Y = model.predict({'x': X, 'f': ['C']*len(X)})
plt.plot(X, Y, 'r-')
Y = model.predict({'x': X, 'f': ['T']*len(X)})
plt.plot(X, Y, 'b-')
In [75]:
draw_scatter(df)
draw_estimate_x_f(x_f_related_model)
plt.xlim(7, 13)
plt.title(u'体サイズ$xi$と施肥効果$fi$を組みこんだ統計モデルによる推定')
Out[75]:
In [44]:
glm('y ~ x + f', data=df, family=sm.families.Poisson(sm.families.links.identity)).fit().summary()
Out[44]:
