この記事では正規分布のパラメータの推定をします．その際に statsmodels と pymc を使用した2通りの方法で行いたいと思います．
学生がメモ代わりに書くものですので，なにか間違いがあればご指摘いただけると幸いです．
データセットの作成
とりあえず今回は簡単にノイズがないデータセットを用意します（正規分布のノイズを載せると全体がことなる正規分布にしたがってしまうため）．
サンプル数が増えることによる流れも見たいのでサンプル数を100, 1000, 10000とします．
import numpy as np
MU = 10.0
SIGMA = 1.0
y0 = np.random.normal(MU, SIGMA, size=100)
y1 = np.random.normal(MU, SIGMA, size=1000)
y2 = np.random.normal(MU, SIGMA, size=10000)
y = [y1, y2, y3]
得られたデータをプロットすると以下のような感じです．
y2になるにつれて正規分布に近づいてきますね．
statsmodelsによる推定
それではstatsmodelsを使用した推定をしていきます．
import statsmodels.api as sm
for i in range(3):
size = 10**(i + 2)
x = [1 for j in range(size)]
glm = sm.GLM(np.array(y[i]), np.array(x), family=sm.families.Gaussian())
res = glm.fit()
print(res.summary())
結果は以下のようになりました．
| サンプル番号 | 平均 | 分散 | 95%信頼区間 (下限) | 95%信頼区間 (上限) |
|---|---|---|---|---|
| 0 | 9.936 | 0.874 | 9.752 | 10.119 |
| 1 | 9.981 | 1.010 | 9.919 | 10.043 |
| 2 | 9.995 | 1.016 | 9.976 | 10.015 |
サンプル数が増えるにつれてちゃんと95%信頼区間が狭まっています（ノイズを乗せていないのであたりまえ）．
Pymcによる推定
pymc3
いよいよPymcで推定したいと思います．
コードは先程と変わって一つのサンプルに対して実行するようにしているので注意してください．
また，pymcはtauを入力とするので分散の値ではないことにも注意してください．
import statsmodels.api as sm
from scipy import optimize
import pymc3 as pm3
with pm3.Model() as model:
# mu, tauの事前分布は一様分布
mu = pm3.Uniform("mu", upper= 10**2, lower= -(10**2))
tau = pm3.Uniform("tau", upper= 10**2, lower= 0)
Y_obs = pm3.Normal("Y_obs", mu=mu, sd=tau, observed=np.array(y3))
with model:
# MAP推定量の計算
start = pm3.find_MAP(fmin=optimize.fmin_powell)
# メトロポリス法のサンプリング
# NUTS, Metropolis, Slice, HamiltonianMCなどがある
step = pm3.Metropolis()
# 5000回サンプリング
trace = pm3.sample(5000, start=start, step=step)
# 2000個を捨てる
pm3.traceplot(trace[2000:])
plt.show()
pm3.summary(trace[2000:])
| サンプル番号 | mu | tau | 95%信用区間 (下限，mu, tauの順) | 95%信用区間 (上限，mu, tauの順) |
|---|---|---|---|---|
| 0 | 9.931 | 0.948 | 9.742, 0.811 | 10.091, 1.080 |
| 1 | 9.979 | 1.007 | 9.915, 0.966 | 10.040, 1.056 |
| 2 | 9.995 | 1.008 | 9.975, 0.993 | 10.016, 1.022 |
また，y2に対するサンプリング結果を載せておきます．
pymc2
pymc2のコードものせておきます．
# モデルの作成
mu = pm.Uniform("mu", upper= 10**2, lower= -(10**2))
tau = pm.Uniform("tau", upper= 10**2, lower= 0)
obs = pm.Normal("obs", mu=mu, tau=tau, value=y3, observed=True)
model = pm.Model([obs, mu, tau])
# MAP推定量の計算
map_ = pm.MAP(model)
map_.fit()
# MCMCの計算
mcmc = pm.MCMC(model)
mcmc.sample(5000, 2000)
# 可視化
pm.Matplot.plot(mcmc.trace("mu"), common_scale=False)
pm.Matplot.plot(mcmc.trace("tau"), common_scale=False)
plt.show()
おわりに
頻度論的な推定とベイズ主義的な推定の手法を試してみました．
statsmodelsとpymcを両方扱っているページがあまりなさそうだったので，自分用のメモをかねて書いてみました．MCMCの各サンプリング手法の中身についてまだいまいち理解できていないので，今後はそこらへんを勉強できたらなと思っています．