Stanのgaussian_dlm_obs()によるトレンドモデル

以前うまくいかなかったStanのgaussian_dlm_obs()によるトレンドモデルだが、モデルを見直してできるようになった。

ナイル川のデータをつかった。Stanによるモデルは以下のように。

data {
  int<lower=0>  N;
  matrix[1, N]  y;
}
transformed data {
  matrix[2, 1]  F;
  matrix[2, 2]  G;
  vector[2]     m0;
  cov_matrix[2] C0;

  F[1, 1] <- 1;
  F[2, 1] <- 0;
  G[1, 1] <- 1;
  G[1, 2] <- 1;
  G[2, 1] <- 0;
  G[2, 2] <- 1;
  m0[1] <- 0;
  m0[2] <- 0;
  C0[1, 1] <- 1.0e+7;
  C0[1, 2] <- 0;
  C0[2, 1] <- 0;
  C0[2, 2] <- 1.0e+7;
}
parameters {
  real<lower=0> sigma[3];
}
transformed parameters {
  vector[1]     V;
  cov_matrix[2] W;

  V[1] <- sigma[1] * sigma[1];
  W[1, 1] <- sigma[2] * sigma[2];
  W[1, 2] <- 0;
  W[2, 1] <- 0;
  W[2, 2] <- sigma[3] * sigma[3];
}
model {
  y ~ gaussian_dlm_obs(F, G, V, W, m0, C0);
  sigma ~ uniform(0, 1.0e+6);
}

前はmatrix[1, 2] F;でエラーになっていたが、matrix[2, 1] F;でうごくようになった。数式は、

というようになるとおもうので、前の方が式にあっていそうなのだが。
【追記】マニュアルをちゃんと読むと、ということで、Fは転置されるのでありました。

Rコード

library(rstan)
library(parallel)

fit.stan2 <- stan(file = "Nile2.stan",
                  data = list(N = length(Nile),
                              y = matrix(c(Nile), nrow = 1)),
                  pars = c("V", "W"),
                  chains = 4,
                  iter = 2000, warmup = 1000, thin = 1)
print(fit.stan2)

V <- get_posterior_mean(fit.stan2, par = "V")[, "mean-all chains"]
W <- get_posterior_mean(fit.stan2, par = "W")[, "mean-all chains"]

library(dlm)

BuildNile <- function(par) {
  dlmModPoly(order = 2, dV = par[1], dW = par[2:3])
}
smoothNile <- dlmSmooth(Nile, BuildNile(c(V, W[1], W[4])))

par(mfrow = c(3, 1))
plot(Nile, las = 1)
plot(dropFirst(smoothNile$s[, 1]), las = 1, ylab = "theta[1]")
plot(dropFirst(smoothNile$s[, 2]), las = 1, ylab = "theta[2]")

結果のグラフ

タグ：STAn