ちょっとStan一辺倒でやってるのも随分効率が悪いなぁと思い始めてきたところに、大仏のオッサンがこんなナイスな記事をupしていたのに今頃気付いたのでした（オッサン気付くの遅くてごめん）。
で、さらに言うと@ibaibabaibai先生からこんなツッコミを頂戴いたしまして。。。
トレンド２階差分で全部正規ってモデルですよね。これでSTAN? 普通はまずカルマンフィルタでは。最小２乗法＋周辺尤度最大化で（超）パラメータ決定でも長さ１００Ｘ３本ならできるのでは。俺なんかまた勘違いしてるのかなあ。http://t.co/RgQLiXY9yv
— baibai (@ibaibabaibai) 2014, 6月 24ということで、一旦Stanシリーズ記事は中休みにして（今後こちらが進展してきたらStan側の記事も書くかも）、ちょっと{dlm}パッケージを触りつつ動的線形モデルについてうだうだ書きながら勉強していこうかと思います。で、今回のシリーズで参考にする資料は以下の2つ。
ベイジアン界隈では高名なHiroki Itô先生のslideshareと、しばらくずーーーっと僕のAmazonのお薦めトップに出てたこの一冊。
- 作者: G.ペトリス,S.ペトローネ,P.カンパニョーリ,和合肇,萩原淳一郎
- 出版社/メーカー: 朝倉書店
- 発売日: 2013/05/08
- メディア: 単行本
- この商品を含むブログを見る
最近カード使用額が多いのでどうしようかなーと迷ってたんですが、結局ポチっちゃいました。ということでこの本を読みながら、そして適宜例題を解いたりしながら、僕自身が学んでいった過程を書いてみようかと思います。
そもそもなぜ{dlm}で「動的」モデルをやるのか
既に「Stanで学ぶ」シリーズ記事で結構複雑なモデルを扱うようになってますが、まぁbaibai先生からもツッコまれた通りで「そもそもStan使うとかやり過ぎなんじゃね」感もあり。
うーむ、たしかにSTANは線形正規モデルに使っても、アルゴリズムにヘッセ行列入ってるから損になる度合いは小さそう。実はそれがメリットだったり。なんか機関銃で人の頭殴ってるみたいな気もするけど。
— baibai (@ibaibabaibai) 2014, 6月 24どうあがいてもやっぱりStanは計算負荷が高いので（笑）、ならもうちょっと計算負荷も低くて汎用パッケージとして{dlm}も用意されている動的線形モデルにトライしてみても良いのかなぁ、と思ったのでした。いや、順番逆だろとツッコまれそうですが。。。
後は、Stanでやろうとすると事後分布が安定しないなぁと思うケースがあったという点。以前階層ベイズで用いたサンプルデータをdatという名前で読み込んで、以下のように加工して目的変数y2を作ってみます。
for (i in 2:100) { a[i]<-0.9*a[i-1]+rnorm(n=1,mean=0,sd=0.00005) b[i]<-0.8*b[i-1]+rnorm(n=1,mean=0,sd=0.0001) c[i]<-0.7*c[i-1]+rnorm(n=1,mean=0,sd=1e-05) d[i]<-0.9*d[i-1]+rnorm(n=1,mean=0,sd=100) } y2<-dat$x1*a+dat$x2*b+dat$x3*c+d+cumsum(c(rep(10,30),rep(20,30),rep(50,40)))
で、こいつを頑張って以下のStanコードで推定してみます。
data {
int<lower=0> N;
real x1[N];
real x2[N];
real x3[N];
real y[N];
}
parameters {
real a[N];
real b[N];
real c[N];
real d[N];
real trend[N];
real s_trend;
real s_q;
real s_a;
real s_b;
real s_c;
real s_d;
real e_a;
real e_b;
real e_c;
real e_d;
}
model {
real q[N];
real cum_trend[N];
trend~normal(30,10);
for (i in 3:N)
trend[i]~normal(2*trend[i-1]-trend[i-2],s_trend);
cum_trend[1]<-trend[1];
for (i in 2:N)
cum_trend[i]<-cum_trend[i-1]+trend[i];
for (i in 1:N)
q[i]<-y[i]-cum_trend[i];
for (i in 2:N) {
a[i]~normal(e_a*a[i-1],s_a);
b[i]~normal(e_b*b[i-1],s_b);
c[i]~normal(e_c*c[i-1],s_c);
d[i]~normal(e_d*d[i-1],s_b);
}
for (i in 1:N)
q[i]~normal(a[i]*x1[i]+b[i]*x2[i]+c[i]*x3[i]+d[i],s_q);
}
出てきた結果がこんな感じ。
Inference for Stan model: hb_trend_dlm. 4 chains, each with iter=1000; warmup=500; thin=1; post-warmup draws per chain=500, total post-warmup draws=2000. mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat e_a -0.06 0.47 0.67 -1.16 -0.41 0.08 0.53 0.68 2 12.50 e_b -0.89 0.40 0.56 -1.67 -1.29 -0.86 -0.46 -0.09 2 13.28 e_c 0.17 0.07 0.11 -0.03 0.04 0.20 0.25 0.30 2 3.52 e_d 0.53 0.59 0.83 -0.61 -0.26 0.65 1.31 1.56 2 17.30 lp__ -1302.82 35.20 54.18 -1432.96 -1332.04 -1293.51 -1260.68 -1231.45 2 4.85 Samples were drawn using NUTS(diag_e) at Fri Jul 25 12:51:21 2014. For each parameter, n_eff is a crude measure of effective sample size, and Rhat is the potential scale reduction factor on split chains (at convergence, Rhat=1).
全然合ってないし、収束もしてないしorz こういう、偏回帰係数自体がそもそもAR過程みたいな、文字通り「動的」にパラメータが変化していくようなケースが鬼門なんですね。こうなるとそもそもいちいちモデル式立ててStanコード書いて、Stan上で毎回C++コンパイラ走らせて、収束しないand/or推定結果が妥当でないと思ったらやり直す、というのを繰り返していたらあっという間にジジイになってしまうので*1、もうちょっと筋の良い方法はないかなぁと思うわけです。この辺の手間もできれば{dlm}で吸収したいという思惑もあったりするのでした。
計量時系列分析ではなく、動的線形モデルでやる理由
このブログでは以前計量時系列分析を徹底して特集したことがありますが、新たに同様の枠組みを持ちながらも違うフレームワークである動的線形モデルに取り組むのには、一応理由があります。
それは単純で、『Rによるベイジアン動的線型モデル』p.78にも書かれていますが「制御工学のアナロジーとして介入操作を用いてシステムの状態を制御する」という観点から、データ分析することが可能だからです*2。特に、状態空間モデルを含むDLM全般は、マーケティングにおいては外生変数を含んだ何かしらのKPI / KGIを説明するモデルとして機能し得るわけで、そこで適切な介入操作をモデルに従って行うことで望むKGI / KPIが得られるのであれば、アドホック分析を担う身としてはこんなに有難い話はないです。
後は、システムの状態を推定するということで欠損値・欠測値に対してもロバストで、さらに言えばオンライン推定も可能*3ということで、アルゴリズム実装方面への道も開けている。この辺を意識しながら勉強していきたいというのが今回のモチベーションです。
最後に
ということで、明らかにStanのようなベイジアンに手を出す前にやっておくべきだったDLMを今さら勉強しようというダメダメな企画なので、詳しい方にはぜひマサカリをどんどん投げていただき、炎上ラーニングを成功させていただきたく。。。orz