時変係数動的線形モデル続き：時変・時不変・OLS線形回帰で比較してみる - 渋谷駅前で働くデータサイエンティストのブログ

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')