六本木で働くデータサイエンティストのブログ

元祖「銀座で働くデータサイエンティスト」です / 道玄坂→銀座→東京→六本木

時変係数動的線形モデルをStanで推定してみる(追記あり)

これはただの備忘録です。目新しい内容は特に何もありません。きちんとした内容を学びたいという方は、先日著者の萩原さんからご恵贈いただいたこちらの書籍で学ばれることをお薦めいたします。MCMCに留まらず、粒子フィルタの実装&実践までカバーしていて素晴らしい教科書だと思います。

で、今回取り上げるのは時変係数から成る動的線形モデルです。マーケティング分析では割と時不変係数モデルが暗黙のうちに採用されることが多いのですが、一方で「おいこれどう見ても特定の特徴量の係数は時変やろ」みたいなケースも時々見かけることがあってやろうとは思っていたものの、単にStanでどう書くか考えるのが面倒で放っていたという(笑)。ということで、今回は真面目にやってみようと思います。


シミュレーションデータを作る


適当に乱数シードを選んで作っておきます。ただし時変係数のところとかはシードをきちんと指定できないので同じようにはいかないと思いますので、皆さんのお手元で適宜調整してください。

> set.seed(1001)
> x1 <- runif(100, 3, 5)
> set.seed(1002)
> x2 <- runif(100, 4, 7)
> set.seed(1003)
> x3 <- runif(100, 8, 15)

> sd1 <- 0.1
> sd2 <- 0.2
> sd3 <- 0.05
> sd_y <- 5

> b1 <- rep(0, 100)
> b1[1] <- 2
> for (i in 2:100){
+     b1[i] <- b1[i-1] + rnorm(1, 0, sd1)
+ }
> b2 <- rep(0, 100)
> b2[1] <- 3
> for (i in 2:100){
+     b2[i] <- b2[i-1] + rnorm(1, 0, sd2)
+ }
> b3 <- rep(0, 100)
> b3[1] <- 0.5
> for (i in 2:100){
+     b3[i] <- b3[i-1] + rnorm(1, 0, sd3)
+ }
> b0 <- 10

> y <- rep(0, 100)
> for (i in 1:100){
+     y[i] <- b0 + b1[i] * x1[i] + b2[i] * x2[i] + b3[i] * x3[i]
+ }
> y <- y + rnorm(100, 0, sd_y)

見ての通りで、基本的には3説明変数と1バイアスから成る動的線形モデルです。ただし個々の説明変数にかかる係数は単位根過程(=ローカル線形トレンド)に従う時変パラメータです。目的変数をプロットするとこうなります。

f:id:TJO:20180519152944p:plain

何かどう見てもローカル線形トレンドか二回差分トレンドのありそうなデータになってしまいましたが、一旦これで良しとします。


Stanで推定する


あとはベタッとStanで書くだけです。本来ならきちんと個々のサンプル間の関係をDAGにして数式で書いておくべきなんでしょうが、ただの備忘録なのでここでは割愛します。Stanコードは以下の通り。例えば"dlm_tv.stan"みたいな名前で保存しておきます。

data {
	int<lower=0> N;
	int<lower=0> M;
	matrix[N, M] X;
	real<lower=0> y[N];
}

parameters {
	real s_q;
	vector<lower=0>[M] s_b;
	matrix<lower=0>[N, M] beta;
	real d_int;
}

model {
	for (i in 2:N){
		for (j in 1:M){
			beta[i, j] ~ normal(beta[i-1, j], s_b[j]);
		}
	}
	for (i in 1:N)
		y[i]~normal(dot_product(X[i,], beta[i,]) + d_int, s_q);
}

これを以下のRコードでkickします。

> mx <- cbind(x1, x2, x3)
> dat <- list(N = nrow(mx), M = ncol(mx), y = y, X = mx)
> library(rstan)
> options(mc.cores = parallel::detectCores())
> rstan_options(auto_write = TRUE)
> fit <- stan(file = 'dlm_tv.stan', data = dat, iter = 1000, chains = 4)
## 中略 ##
> slength <- 2000
> fit.smp <- extract(fit)
> tmp <- density(fit.smp$d_int)
> d_int <- tmp$x[tmp$y == max(tmp$y)]
> tmp <- density(fit.smp$s_q)
> s_q <- tmp$x[tmp$y == max(tmp$y)]
> s_beta <- rep(0, ncol(mx))
> for (i in 1:ncol(mx)) {
#     tmp <- density(fit.smp$beta[(slength*(i-1)+1):(slength*i)])
+     tmp <- density(fit.smp$s_b[(slength*(i-1)+1):(slength*i)]) # 萩原さんからバグをご指摘いただいたのでこちらに修正
+     s_beta[i] <- tmp$x[tmp$y == max(tmp$y)]
+ }
> tmp <- fit.smp$beta
> beta <- matrix(0, nrow = nrow(mx), ncol = ncol(mx))
> for (i in 1:nrow(mx)){
+     for (j in 1:ncol(mx)){
+         tmpb <- density(tmp[, i, j])
+         beta[i, j] <- tmpb$x[tmpb$y == max(tmpb$y)]
+     }
+ }
> y_pred <- rep(0, 100)
> for (i in 1:100){
+     y_pred[i] <- d_int + beta[i, 1] * x1[i] + beta[i, 2] * x2[i] + beta[i, 3] * x3[i]
+ }
> matplot(cbind(y, y_pred), type = 'l', xlab = '', ylab = '', lwd = 2)
> legend('topright', legend = c('Data', 'Fitted'), col = c(1,2), lty = c(1,2), lwd = 2)

f:id:TJO:20180519153701p:plain

見た感じフィッティング感は悪くないようです。では時変係数のフィッティングを見てみましょう。

> par(mfrow = c(3, 1))
> matplot(cbind(b1, beta[, 1]), type = 'l', xlab = '', ylab = '', lwd = 2, main = 'Beta 1')
> matplot(cbind(b2, beta[, 2]), type = 'l', xlab = '', ylab = '', lwd = 2, main = 'Beta 2')
> matplot(cbind(b3, beta[, 3]), type = 'l', xlab = '', ylab = '', lwd = 2, main = 'Beta 3')

f:id:TJO:20180519154211p:plain

これ全然合ってない気がしますねorz まぁ"All models are wrong, but some are useful"ということで時変係数モデルのようなフリーパラメータの多いモデルだとこういうこともあるのかなと。最後にMCMCサンプリングの収束具合を見ておきます。

> options(max.print = 10000)
> fit
Inference for Stan model: dlm_tv.
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
s_q           4.22    0.02  0.41   3.49   3.93   4.21   4.50   5.08   360 1.01
s_b[1]        0.18    0.03  0.12   0.02   0.08   0.16   0.25   0.46    18 1.23
s_b[2]        0.14    0.03  0.09   0.02   0.07   0.13   0.20   0.35     7 1.30
s_b[3]        0.08    0.01  0.04   0.02   0.06   0.08   0.11   0.16    30 1.09
beta[1,1]     2.90    0.16  1.04   1.13   2.16   2.82   3.53   5.22    40 1.06
beta[1,2]     2.48    0.09  0.85   0.82   1.90   2.46   3.06   4.09    84 1.05
beta[1,3]     1.09    0.03  0.38   0.42   0.82   1.08   1.35   1.85   135 1.03
beta[2,1]     2.89    0.17  1.04   1.13   2.14   2.80   3.49   5.23    38 1.06
beta[2,2]     2.46    0.09  0.83   0.83   1.91   2.44   3.03   4.05    86 1.05
beta[2,3]     1.08    0.03  0.37   0.43   0.81   1.07   1.32   1.82   129 1.02
## 中略 ##
beta[99,1]    2.44    0.06  0.89   0.89   1.81   2.37   3.00   4.41   251 1.00
beta[99,2]    1.82    0.06  0.69   0.53   1.36   1.80   2.29   3.21   135 1.02
beta[99,3]    0.56    0.02  0.29   0.10   0.34   0.53   0.76   1.17   151 1.01
beta[100,1]   2.42    0.06  0.92   0.82   1.76   2.36   3.00   4.47   244 1.00
beta[100,2]   1.80    0.06  0.70   0.46   1.32   1.78   2.29   3.24   135 1.02
beta[100,3]   0.56    0.02  0.29   0.09   0.34   0.53   0.75   1.17   144 1.01
d_int         4.29    0.45  4.86  -5.73   1.07   4.49   7.64  13.44   115 1.02
lp__        442.13   18.16 84.74 308.06 379.70 430.10 496.34 631.20    22 1.11

Samples were drawn using NUTS(diag_e) at Sat May 19 15:16:18 2018.
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).

Rhatは大体1.00に近いあたりをどのパラメータでも推移していてそこまで悪くはなさそうです。最後に{coda}でMCMCサンプルのヒストグラムを見てみましょう。

> library(coda)
> fit.coda<-mcmc.list(lapply(1:ncol(fit),function(x) mcmc(as.array(fit)[,x,])))
> plot(fit.coda[,5:7])
> plot(fit.coda[,302:304])
> plot(fit.coda[,305:306])

f:id:TJO:20180519155110p:plain
f:id:TJO:20180519155131p:plain
f:id:TJO:20180519155149p:plain

ヒストグラムを見た感じもそこまで収束具合は悪くなさそうです。とりあえず、今回のやり方で一応時変係数動的線形モデルの推定がStanで出来るということは分かりました。後はきちんとこれが適用できて尚且つフィッティングも良くてさらに汎化性能も確保できるようなやり方と、データの組み合わせを探す感じですかね。。。


追記


『基礎からわかる時系列分析』の著者、萩原さんから以下のように{dlm}パッケージを使った場合のフォローアップ記事を頂戴しました。萩原さん、有難うございます。この程度の計算なら&特に事後分布が欲しいわけではないのであれば、{dlm}パッケージを使った方が圧倒的に速いです(試してみてやっぱりなぁ…となりました)。

また上のRコード内のバグもご指摘いただきましたので、1点修正してあります。大変失礼いたしましたm(_ _)m


ところで萩原さんからは、

この結果を見ると、モデルのパラメータを真値に設定すると状態(平均値)の推定精度は若干向上するものの、限界はある印象を受けます。 従って、回帰係数の推定精度を追求する観点では、今回のデータは元々やや手ごわいものになっているのではないかと感じます。

という感想を頂戴しているのですが、後になって見返してみたら「そりゃそうだよな」と思ってしまいました。理由は簡単で、3つの真の時変係数同士を比べてみると

> matplot(scale(cbind(b1, b2, b3)), type = 'l', lty = 1, ylab = '')
> cor(cbind(b1,b2,b3))
           b1         b2         b3
b1  1.0000000 -0.6065728 -0.4004846
b2 -0.6065728  1.0000000  0.8616015
b3 -0.4004846  0.8616015  1.0000000

f:id:TJO:20180601121714p:plain
見ての通りb2とb3とで相関が高く、多重共線性が高くなるのと同様の現象を起こしていたからですorz もうちょっと適切なシミュレーションデータを作るべきだったと反省しております。。。