DLMをやるとどうもオーバーフィッティングしているように見える気がしたので 自分用のメモ＋可視化． いろいろ考えたあと実際の数字をいろいろ見てるとシステムノイズを工夫すると 予測との乖離から外的要因が見えるようになって視界がクリアになった気がした． クリエイティブの変化とか，消費税のあれとか． どんなに頑張ってもすべての外的要因をモデリングすることなんてできなそうなので それを定数項が吸収するくらいなら思い切って予測を外してほしいし ノイズみたいなのを拾いたくもない

require(rstan)
require(dplyr)
require(magrittr)
require(ggvis)
require(doParallel)

データの作成

本当は神様しか知らないデータ達

T <- 100
a0 <- rnorm(T,0,0.1) %>% cumsum()
a1 <- rnorm(T,0,0.15) %>% cumsum()
a2 <- rnorm(T,0,0.02) %>% cumsum()

x1 <- rpois(T,20)
x2 <- rpois(T,20)

y <- a0 + a1 * x1 + a2 * x2



dat.list <- list(x1=x1, x2=x2, y=y, T=T )

それぞれのStan

stancode <- '
data{
  int T;
  vector[T] x1;
  vector[T] x2;
  vector[T] y;
}
parameters{
  vector[T] a0;
  vector[T] a1;
  vector[T] a2;

  real<lower=0> obs_y;
  real<lower=0> sn_a0;
  real<lower=0> sn_a1;
  real<lower=0> sn_a2;
  
}
transformed parameters{
  vector[T] yp;
  yp <- a0 + a1 .* x1 + a2 .* x2;
}
model{
# priors
  a0[1] ~ normal(0,100);
  a1[1] ~ normal(0,100);
  a2[1] ~ normal(0,100);

  obs_y ~ normal(0,100);
  sn_a0 ~ normal(0,100);
  sn_a1 ~ normal(0,100);
  sn_a2 ~ normal(0,100);

  for(t in 2:T){
    a0[t] ~ normal(a0[t-1],sn_a0);
    a1[t] ~ normal(a1[t-1],sn_a1);
    a2[t] ~ normal(a2[t-1],sn_a2);

  }

# likelihood
  y ~ normal(yp,obs_y);

}

'

stancode2 <- '
data{
  int T;
  vector[T] x1;
  vector[T] x2;
  vector[T] y;
}
parameters{
  vector[T] a0;
  vector[T] a1;
  vector[T] a2;

  real<lower=0> obs_y;
  real<lower=0> sn_a0;
  real<lower=0> sn_a1;
  real<lower=0> sn_a2;
  
}
transformed parameters{
  vector[T] yp;
  yp <- a0 + a1 .* x1 + a2 .* x2;
}
model{
# priors
  a0[1] ~ normal(0,100);
  a1[1] ~ normal(0,100);
  a2[1] ~ normal(0,100);

  a0[2] - a0[1] ~ normal(0,100);
  a1[2] - a1[1] ~ normal(0,100);
  a2[2] - a2[1] ~ normal(0,100);

  obs_y ~ normal(0,100);
  sn_a0 ~ normal(0,100);
  sn_a1 ~ normal(0,100);
  sn_a2 ~ normal(0,100);

  for(t in 3:T){
    a0[t] - a0[t-1] ~ normal(a0[t-1] - a0[t-2],sn_a0);
    a1[t] - a1[t-1] ~ normal(a1[t-1] - a1[t-2],sn_a1);
    a2[t] - a2[t-1] ~ normal(a2[t-1] - a2[t-2],sn_a2);

  }

# likelihood
  y ~ normal(yp,obs_y);


}

'

stanによるシミュレーション

chains <- 2
cl <- makeCluster( chains)
registerDoParallel( cl) #並列処理開始の呪文
sflist <- foreach( i=1:chains) %dopar% {
  require( rstan)
  stan(model_code=stancode, data=dat.list, chains=1, iter=500)
}

fit <- sflist2stanfit( sflist)
plot(fit)

a0の95%信用区間と実績

a0p <- rstan::extract(fit)$a0 %>%
  apply(2,median)

ci_a0p <- rstan::extract(fit)$a0 %>%
  apply(2,function(x){quantile(x,probs = c(0.025,0.975))})

a0.df <- data.frame(a0=a0, a0p=a0p, a0pu=ci_a0p[2,], a0pl=ci_a0p[1,],x=1:T)

a0.df %>%
  ggvis( opacity:= 1) %>%
    layer_ribbons(~x,~a0pu,y2 = ~a0pl, opacity := 0.3 ,fill := 'green') %>%
  layer_paths(~x,~a0,stroke := 'red') %>%
  layer_paths(~x,~a0p,stroke := 'green')

a1の95%信用区間と実績

a1p <- rstan::extract(fit)$a1 %>%
  apply(2,median)

ci_a1p <- rstan::extract(fit)$a1 %>%
  apply(2,function(x){quantile(x,probs = c(0.025,0.975))})

a1.df <- data.frame(a1=a1, a1p=a1p, a1pu=ci_a1p[2,], a1pl=ci_a1p[1,],x=1:T)

a1.df %>%
  ggvis( opacity:= 1) %>%
    layer_ribbons(~x,~a1pu,y2 = ~a1pl, opacity := 0.3 ,fill := 'green') %>%
  layer_paths(~x,~a1,stroke := 'red') %>%
  layer_paths(~x,~a1p,stroke := 'green')

a2の95%信用区間と実績

a2p <- rstan::extract(fit)$a2 %>%
  apply(2,median)

ci_a2p <- rstan::extract(fit)$a2 %>%
  apply(2,function(x){quantile(x,probs = c(0.025,0.975))})

a2.df <- data.frame(a2=a2, a2p=a2p, a2pu=ci_a2p[2,], a2pl=ci_a2p[1,],x=1:T)

a2.df %>%
  ggvis( opacity:= 1) %>%
    layer_ribbons(~x,~a2pu,y2 = ~a2pl, opacity := 0.3 ,fill := 'green') %>%
  layer_paths(~x,~a2,stroke := 'red') %>%
  layer_paths(~x,~a2p,stroke := 'green')

yの95%信用区間と実績

yp <- rstan::extract(fit)$yp %>%
  apply(2,median)

ci_yp <- rstan::extract(fit)$yp %>%
  apply(2,function(x){quantile(x,probs = c(0.025,0.975))})

y.df <- data.frame(y=y, yp=yp, ypu=ci_yp[2,], ypl=ci_yp[1,],x=1:T)

y.df %>%
  ggvis( opacity:= 1) %>%
    layer_ribbons(~x,~ypu,y2 = ~ypl, opacity := 0.3 ,fill := 'green') %>%
  layer_paths(~x,~y,stroke := 'red') %>%
  layer_paths(~x,~yp,stroke := 'green')

chains <- 2
cl <- makeCluster( chains)
registerDoParallel( cl) #並列処理開始の呪文
sflist <- foreach( i=1:chains) %dopar% {
  require( rstan)
  stan(model_code=stancode2, data=dat.list, chains=1, iter=500)
}
## Warning: closing unused connection 6 (<-localhost:11408)
## Warning: closing unused connection 5 (<-localhost:11408)
fit2 <- sflist2stanfit( sflist)
plot(fit2)

stopCluster(cl)

a0の95%信用区間と実績(2回階差)

a0p <- rstan::extract(fit2)$a0 %>%
  apply(2,median)

ci_a0p <- rstan::extract(fit2)$a0 %>%
  apply(2,function(x){quantile(x,probs = c(0.025,0.975))})

a0.df <- data.frame(a0=a0, a0p=a0p, a0pu=ci_a0p[2,], a0pl=ci_a0p[1,],x=1:T)

a0.df %>%
  ggvis( opacity:= 1) %>%
    layer_ribbons(~x,~a0pu,y2 = ~a0pl, opacity := 0.3 ,fill := 'green') %>%
  layer_paths(~x,~a0,stroke := 'red') %>%
  layer_paths(~x,~a0p,stroke := 'green')

a1の95%信用区間と実績(2回階差)

a1p <- rstan::extract(fit2)$a1 %>%
  apply(2,median)

ci_a1p <- rstan::extract(fit2)$a1 %>%
  apply(2,function(x){quantile(x,probs = c(0.025,0.975))})

a1.df <- data.frame(a1=a1, a1p=a1p, a1pu=ci_a1p[2,], a1pl=ci_a1p[1,],x=1:T)

a1.df %>%
  ggvis( opacity:= 1) %>%
    layer_ribbons(~x,~a1pu,y2 = ~a1pl, opacity := 0.3 ,fill := 'green') %>%
  layer_paths(~x,~a1,stroke := 'red') %>%
  layer_paths(~x,~a1p,stroke := 'green')

a2の95%信用区間と実績(2回階差)

a2p <- rstan::extract(fit2)$a2 %>%
  apply(2,median)

ci_a2p <- rstan::extract(fit2)$a2 %>%
  apply(2,function(x){quantile(x,probs = c(0.025,0.975))})

a2.df <- data.frame(a2=a2, a2p=a2p, a2pu=ci_a2p[2,], a2pl=ci_a2p[1,],x=1:T)

a2.df %>%
  ggvis( opacity:= 1) %>%
    layer_ribbons(~x,~a2pu,y2 = ~a2pl, opacity := 0.3 ,fill := 'green') %>%
  layer_paths(~x,~a2,stroke := 'red') %>%
  layer_paths(~x,~a2p,stroke := 'green')

yの95%信用区間と実績(2回階差)

yp <- rstan::extract(fit2)$yp %>%
  apply(2,median)

ci_yp <- rstan::extract(fit2)$yp %>%
  apply(2,function(x){quantile(x,probs = c(0.025,0.975))})

y.df <- data.frame(y=y, yp=yp, ypu=ci_yp[2,], ypl=ci_yp[1,],x=1:T)

y.df %>%
  ggvis( opacity:= 1) %>%
    layer_ribbons(~x,~ypu,y2 = ~ypl, opacity := 0.3 ,fill := 'green') %>%
  layer_paths(~x,~y,stroke := 'red') %>%
  layer_paths(~x,~yp,stroke := 'green')