廿TT

譬如水怙牛過窓櫺 頭角四蹄都過了 因甚麼尾巴過不得

pymc3 で状態空間(ローカルレベル)モデル

Python ハローワールドです。

Rodeo(Yhat End-to-End Data Science Platform: Rodeo)っていう RStudio みたいなやつをいれました。

なかなか快適です。

import numpy as np
import pymc3 as pm
import matplotlib.pyplot as plt
import statistics as stat
import pickle as pkl
from scipy import optimize

np.random.seed(12)
y_obs = np.cumsum(np.random.normal(size=100))

plt.plot(y_obs)

とりあえずランダムウォークで適当な系列をつくりました。

f:id:abrahamcow:20170329222312p:plain

ローカルレベルモデルを pymc3 で書くとこんな感じになりました。

with pm.Model() as model1:
    tau = pm.Uniform('tau',lower=0,upper=1e+7)
    states = [pm.Normal("s0", mu=0, tau=tau)]
    for i in range(1, 100):
        states.append(pm.Normal(name="s" + str(i), mu=states[i-1], tau=tau))
        
    y = pm.Normal("y", mu=states, tau=tau, observed=y_obs)

状態(state)を一回一回アペンドしていってます。あまり賢いやりかたではないような気もしますが、他の方法が思いつかなかった。

まずは MAP 推定してプロットしてみます。

map_estimate = pm.find_MAP(model=model1)

state_map = [map_estimate["s%i" % i] for i in range(0,100)]

with plt.style.context('fivethirtyeight'):
    plt.plot(y_obs, linewidth=2)
    plt.plot(state_map, linewidth=2)

f:id:abrahamcow:20170329222836p:plain

はい。

続いて MCMC です。バーンインは問答無用で二万回になるみたいです。2000は本番(バーンイン期間の後)の回数です。

with model1:
    trace1 = pm.sample(2000,start=map_estimate,chain=3,njobs=3)

しばらく待つと終わります。

トレースプロットはこんな風に確認できます。

pm.traceplot(trace1,varnames=["s0","tau"])

f:id:abrahamcow:20170329223122p:plain

事後期待値と95%ベイズ信頼区間を for 文で取り出します。

これもあまり賢くない気がします。

state_eap = [stat.mean(trace1["s%i" % i]) for i in range(0,100)]
upper = [np.percentile(trace1["s%i" % i],2.5) for i in range(0,100)]
lower = [np.percentile(trace1["s%i" % i],97.5) for i in range(0,100)]


with plt.style.context('fivethirtyeight'):
    xseq = range(0,100)
    plt.plot(xseq,y_obs, linewidth=2,color="black")
    plt.plot(xseq,state_eap, linewidth=3,color="black",linestyle=":")
    plt.fill_between(xseq,lower,upper,alpha=0.5,color="grey")

灰色の帯が95%ベイズ信頼区間です。

f:id:abrahamcow:20170329223301p:plain

以上です。

参考文献は岩波データサイエンスの一巻です。

岩波データサイエンス Vol.1

岩波データサイエンス Vol.1

ローカルレベルモデルのなんたるか、pymc のなんたるかはこちらを参照してください。

pymc3 を触ってみた感想

  • 離散パラメータを扱えるという点に魅力を感じた。
  • でも基本的にはこれまで通り R と Stan でよさそう。