2週間前にふと思い立ってこんなことを試してみたわけですが。
よくよく考えてみたら「データを生成した真のモデルが時変係数＆モデル推定も時変係数」でやってみた結果を並べただけで、これを（例えば）時不変係数モデルで推定してしまった場合や単なるOLS線形回帰で推定した場合との比較みたいな、ありがちなケースを試してみるのを忘れていたのでした。
ということで、単にそれらをやってみます。なお今回も面倒なので交差検証は入れていません。興味のある方は以下の過去記事をご参照ください。
Stanコード
まず、前回同様以下の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);
}
次に、以下のStanコードを"dlm_const.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] beta;
real d_int;
}
model {
for (i in 1:N)
y[i]~normal(dot_product(X[i,], beta) + d_int, s_q);
}
下準備はこれでおしまいです。
Rコードでkickする
あとは以下のRコードでまとめてkickするだけです。そう言えばこの記事書いてて気付いたんですが、MCMCサンプル抽出してきてサンプルの最頻値に対応するパラメータ選ぶところってもっと簡単に書けるんですよねorz そしてそういうMAP推定しか使わないのであればrstan::optimizing使った方が早いので、今度からそちらにしようかなと。。。
# データ生成 > 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) # 時変係数モデル > 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) > fit1 <- stan(file = 'dlm_tv.stan', data = dat, iter = 1000, chains = 4) ## 中略 ## > slength <- 2000 > fit1.smp <- extract(fit1) > tmp <- density(fit1.smp$d_int) > d_int1 <- tmp$x[tmp$y == max(tmp$y)] > tmp <- density(fit1.smp$s_q) > s_q1 <- tmp$x[tmp$y == max(tmp$y)] > s_beta1 <- rep(0, ncol(mx)) > for (i in 1:ncol(mx)) { + tmp <- density(fit1.smp$beta[(slength*(i-1)+1):(slength*i)]) + s_beta1[i] <- tmp$x[tmp$y == max(tmp$y)] + } > tmp <- fit1.smp$beta > beta1 <- 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]) + beta1[i, j] <- tmpb$x[tmpb$y == max(tmpb$y)] + } + } > y_pred1 <- rep(0, 100) > for (i in 1:100){ + y_pred1[i] <- d_int1 + beta1[i, 1] * x1[i] + beta1[i, 2] * x2[i] + beta1[i, 3] * x3[i] + } # 時不変係数モデル > fit2 <- stan(file = 'dlm_const.stan', data = dat, iter = 1000, chains = 4) > tmp <- density(fit2.smp$d_int) > d_int2 <- tmp$x[tmp$y == max(tmp$y)] > tmp <- density(fit2.smp$s_q) > s_q2 <- tmp$x[tmp$y == max(tmp$y)] > beta2 <- rep(0, ncol(mx)) > for (i in 1:ncol(mx)) { + tmp <- density(fit2.smp$beta[(slength*(i-1)+1):(slength*i)]) + beta2[i] <- tmp$x[tmp$y == max(tmp$y)] + } > y_pred2 <- rep(0, 100) > for (i in 1:100){ + y_pred2[i] <- d_int2 + beta2[1] * x1[i] + beta2[2] * x2[i] + beta2[3] * x3[i] + } # OLS線形回帰 > d <- as.data.frame(cbind(y, mx)) > names(d) <- c('y', 'ad1', 'ad2', 'ad3') > d.lm <- lm(y~., d) > y_pred3 <- predict(d.lm, newdata = d) > summary(d.lm) Call: lm(formula = y ~ ., data = d) Residuals: Min 1Q Median 3Q Max -25.116 -7.017 2.210 7.238 17.149 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 7.1526 11.7294 0.610 0.5434 ad1 0.7300 1.6445 0.444 0.6581 ad2 2.3240 1.1916 1.950 0.0540 . ad3 0.9748 0.4960 1.965 0.0523 . --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 9.92 on 96 degrees of freedom Multiple R-squared: 0.07447, Adjusted R-squared: 0.04554 F-statistic: 2.575 on 3 and 96 DF, p-value: 0.0584 > par(mfrow=c(1,1)) > matplot(cbind(y, y_pred1, y_pred2, y_pred3), type = 'l', xlab = '', ylab = '', col = c(1,2,3,4), lty = 1, lwd = 2) > legend('bottomleft', legend = c('Data', 'Time-variant', 'Const', 'OLS'), col = c(1,2,3,4), lty = 1, lwd = 2, ncol = 2) > par(mfrow = c(3, 1)) > matplot(cbind(b1, beta1[,1], rep(beta2[1], nrow(mx)), rep(d.lm$coefficients[2], nrow(mx))), type = 'l', xlab = '', ylab = '', lty = 1, lwd = 2, main = 'Beta 1') > matplot(cbind(b2, beta1[,2], rep(beta2[2], nrow(mx)), rep(d.lm$coefficients[3], nrow(mx))), type = 'l', xlab = '', ylab = '', lty = 1, lwd = 2, main = 'Beta 2') > matplot(cbind(b3, beta1[,3], rep(beta2[3], nrow(mx)), rep(d.lm$coefficients[4], nrow(mx))), type = 'l', xlab = '', ylab = '', lty = 1, lwd = 2, main = 'Beta 3')
全体として、時変係数モデルの当てはまりが最も良くて、時不変係数モデルもOLS線形回帰モデルもあまり当てはまりが良くないことが見て取れます。特に真の回帰係数との比較において、時不変係数モデルとOLS線形回帰では真の値からかなり外れた値に推定されてしまったという結果にはちょっと意外な印象を覚えました。
感想
当たり前っちゃ当たり前なんですが、真のモデルが時変係数の場合に時不変係数モデルを選択すると、そもそも学習データへの当てはまり自体が悪くなるということが実際にやってみて良く分かりました。現場からは以上です。