Pythonでポアソン分布回帰モデルまわりの書き方

データ解析のための統計モデリング入門――一般化線形モデル・階層ベイズモデル・MCMC (確率と情報の科学)

pythonで書くのにはどうすれば良いかを確認しながら読んでいる。

ポアソン分布を書く

from scipy.stats import poisson

x = np.arange(0, 20, 1)
for mu in [3.5, 7.7, 15.1]:
    plt.plot(x, poisson.pmf(x, mu), label='{}'.format(mu))
        
plt.legend()   
plt.show()

ポアソン分布に従ったランダムな値を取得

平均 のポアソン分布に従った値を100個取得する。

import numpy as np
np.random.poisson(lam=7, size=100)

確認

import matplotlib.pyplot as plt
plt.hist(data, bins=15)

histgramの値を整形して、グラフを作成する

n, bins = np.histogram(data, bins=15)

plt.scatter(x=bins[:-1], y=n, linewidths=1)
plt.plot(bins[:-1], n, linewidth=1, linestyle="dashed")

最尤推定値

グラフ化して、ざっくり最尤推定値がどこになるかをみる

from scipy import stats

l = np.arange(0.1, 20, 0.1)
# stats.poisson.logpmf(): 対数尤度関数
logL = np.array([sum(stats.poisson.logpmf(data, i)) for i in l])

x = np.linspace(0, 20, len(logL))
plt.plot(x, logL)
plt.axvline(x=x[logL.argmax()], color='red', linestyle='dotted')

に近い値になる。

ポアソン回帰モデル

まずはデータ取得。

生態学データ解析 - 本/データ解析のための統計モデリング入門から。

df = pd.read_csv('http://hosho.ees.hokudai.ac.jp/~kubo/stat/iwanamibook/fig/poisson/data3a.csv')
df.head()

	y	x	f
0	6	8.31	C
1	6	9.44	C
2	6	9.50	C
3	12	9.07	C
4	10	10.16	C

因子型の数値変換

f列は、TかCが埋まっている。

df.f.head()
0    C
1    C
2    C
3    C
4    C
Name: f, dtype: object

Rだと、因子型が文字列でもよろしくやってくれるが、Pythonの場合は変換をする必要がある。 pandas.get_dummies で変換できる

df['f'] = pd.get_dummies(df['f'])['T']
df.f.head()
0    0
1    0
2    0
3    0
4    0
Name: f, dtype: uint8

statsmodels.api.GLM を使う

import statsmodels.api as sm

x = sm.add_constant(df.x)
fit = sm.GLM(df.y, x, family=sm.families.Poisson()).fit()
fit.summary()

こちらは自動で

statsmodels.formula.api.poisson を使う

import statsmodels.formula.api as smf

fit = smf.poisson("y ~ x", data=df).fit()
fit.summary()

線形予測子の書き方

それぞれ、書き方がある

	線形予測子	smf.poisson のformula	statsmodels.api.add_constantのdata
一定		y ~ 1	np.repeat(1, len(df.x))
x		y ~ x	df.x
x + f		y ~ x + f	df[['x', 'f']]

AICの取得

statsmodels.api.GLM も statsmodels.formula.api.poisson も

result.aic
# 474.7725015397216

参考

• Pythonで「データ解析のための統計モデリング入門」：第2章
• データ解析のための統計モデリング入門 一般化線形モデル（GLM）
• データ解析のための統計モデリング入門