Python
statsmodels
pymc

正規分布のパラメータ推定をpymcとstatsmodelsでやってみたメモ

More than 1 year has passed since last update.

この記事では正規分布のパラメータの推定をします.その際に statsmodelspymc を使用した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]

得られたデータをプロットすると以下のような感じです.

fig1.png

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に対するサンプリング結果を載せておきます.

fig2.png

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の各サンプリング手法の中身についてまだいまいち理解できていないので,今後はそこらへんを勉強できたらなと思っています.