前回の記事では盛大にトレンドつきモデルの式をトレンド累積値でモデリングしないという間抜けなことをしてしまい大変失礼しましたorz
さて、階層ベイズモデルでは際限なく色々な要素を足していくことで、果てしなく複雑っぽいモデルを作ることができるわけです。言い換えれば、線形加算可能な要素である限りは、何をどれほど足しても一応推定は可能なはず。@berobero11さんのところでもこんな試みをされていましたね。
ちなみに今現在@berobero11さんは@smrmktさんとwebスクレイピング＋CAR modelを駆使した「割安賃貸物件探し」という大変高度なネタを展開されているようで、未だにこの程度の初歩的な話でゼーハー言ってる僕など全く手の届かないところにいらっしゃるようで。。。僕は細々とやります（泣）。
で、話を戻すとこの手の線形加算系のモデルでお決まりのパターンとして入ってくるのが季節調整です（以前stl関数を用いて季節調整するパターンを記事にしてます）。今回はStanでその辺をモデリングする方法をさっくり書いてみます。
サンプルデータについて
前回と同じようにGitHubからサンプルデータを落としてきて、dとかいう名前でインポートしておきましょう。そしたら、今度は以下のように落としてきたデータに小細工をします。
> wk<-rep(c(1000,700,300,500,600,900,1300),15) > wk<-wk[1:100] > y2<-round(d$y+wk+rnorm(n=100,mean=0,sd=100),0) > dat2<-list(N=100,y=y2,x1=d$x1,x2=d$x2,x3=d$x3)
これで、1週間周期（7点ごと）の季節成分と正規分布ノイズが目的変数yに加わりました。新たな目的変数y2を含むデータセットdat2の内訳はこんな感じです。
パッと見た目には分かりませんが、結構大きめの季節成分を加えこんであります。なので、以前のモデルだと推定精度が悪くなることが予想されます。実際、前回の記事と同じトレンド累積項のみを与えたモデルで当てはめた結果は以下の通り。
相関係数を計算すると0.9291ということで、決して悪くはないんですが今一つという感じがしますね。ということで、実際にStanコードを書いて季節調整をかけてみましょう。
Stanで当てはめてみる
と言っても、やることは簡単です。要は「7日ごとに同じ値がやってきて（例えば）正規分布に従ったバラツキが乗る」という項を作って、ただ加算するだけです。なのでモデルとしては以下のような感じになります。
←ココが今回加わったところ
追記
正しくは以下の通りです。
これを記述するStanコードは以下の通りです。GitHubからソースコードを落とすこともできます。
data {
int<lower=0> N;
real<lower=0> x1[N];
real<lower=0> x2[N];
real<lower=0> x3[N];
real<lower=0> y[N];
}
parameters {
real wk[N];
real trend[N];
real s_trend;
real s_q;
real s_wk;
real<lower=0> a;
real<lower=0> b;
real<lower=0> c;
real d;
}
model {
real q[N];
real cum_trend[N];
trend~normal(30,10);
wk~normal(500,400);
for (i in 8:N)
wk[i]~normal(wk[i-7],s_wk); // 7日周期季節成分
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]-wk[i]-cum_trend[i]; // 要素分解して残りを線形モデルへ
for (i in 1:N)
q[i]~normal(a*x1[i]+b*x2[i]+c*x3[i]+d,s_q);
}
ここまで来ればいかようにもベイジアンモデルを推定できますが、例えばこんな感じのRコードで算出できます。
> require("rstan") > fit_wk<-stan(file='hb_trend_cum_wk.stan',data=dat2,iter=1000,chains=4) > fit_wk.smp<-extract(fit_wk) > dens_wk_a<-density(fit_wk.smp$a) > dens_wk_b<-density(fit_wk.smp$b) > dens_wk_c<-density(fit_wk.smp$c) > dens_wk_d<-density(fit_wk.smp$d) > a_est_wk<-dens_wk_a$x[dens_wk_a$y==max(dens_wk_a$y)] > b_est_wk<-dens_wk_b$x[dens_wk_b$y==max(dens_wk_b$y)] > c_est_wk<-dens_wk_c$x[dens_wk_c$y==max(dens_wk_c$y)] > d_est_wk<-dens_wk_d$x[dens_wk_d$y==max(dens_wk_d$y)] > trend_est_wk<-rep(0,100) > for (i in 1:100) { + tmp<-density(fit_wk.smp$trend[,i]) + trend_est_wk[i]<-tmp$x[tmp$y==max(tmp$y)] + } > week_est_wk<-rep(0,100) > for (i in 1:100) { + tmp<-density(fit_wk.smp$wk[,i]) + week_est_wk[i]<-tmp$x[tmp$y==max(tmp$y)] + } > pred_wk<-a_est_wk*d$x1+b_est_wk*d$x2+c_est_wk*d$x3+d_est_wk+cumsum(trend_est_wk)+week_est_wk > matplot(cbind(y2,pred_wk),type='l',lty=1,lwd=c(3,2),col=c('red','blue'),ylab="CV") > legend(75,4500,c("Data","Predicted"),col=c('red','blue'),lty=c(1,1),lwd=c(3,2),cex=1.2)
当てはめた結果が以下の通り。実用上は文句なしといったところでしょうか。
ちなみに相関を計算すると0.9814ということで、当てはめとしてはむしろオーバーフィッティングが気になる感じですね（汗）。季節成分について、生成された元のデータとStanで推定したデータを比べたものが以下の図です。
大体合ってるというか、ゲタの分だけずれている程度で残りは悪くないと思います。最後に、季節調整ありモデルとなしモデルとで当てはまり具合を比べるとこうなります。
確かに季節成分に綺麗に追従できているか否かで、モデルの当てはまり具合が変わっていることが見て取れると思います。ただし、実は偏回帰係数にこの辺を推定したしわ寄せが来てしまっていて、
> a_est_wk [1] 0.0001466718 > b_est_wk [1] 3.052333e-05 > c_est_wk [1] 4.496809e-05 > d_est_wk [1] 1554.647
となっており、bとdがずれてしまっています。dがずれるのは季節調整のゲタを合わせたことになっているだけなので構わないんですが、bがずれるのは困るなぁ。。。ということで皆様からのツッコミをお待ちしておりますm(_ _)m
*1:ココが前回の記事で@motivic_氏に突っ込まれて直したところ
