状態空間モデルをためす

状態空間モデルの勉強。やはり自分でコードを書いてみないと。

まずはPetris and Petrone (2011)にあった例から。

data(Nile)

## dlm
##
## Petris G., Petrone, S. (2011) State Space Models in R
## Journal of Statistical Software 41(4)
## http://www.jstatsoft.org/v41/i04
##

library(dlm)

buildNile <- function(theta) {
  dlmModPoly(order = 1, dV = theta[1], dW = theta[2])
}
fit <- dlmMLE(Nile, parm = c(100, 2), buildNile, lower = rep(1e-4, 2))
modNile <- buildNile(fit$par)
smoothNile <- dlmSmooth(Nile, modNile)

plot(Nile, type ='o', las = 1)
lines(dropFirst(smoothNile$s), col = "red")

つぎにベイズモデリングでやってみる。自作JAGS呼び出し関数をつかう。このモデルでAR(1)ということになるだろうか。

## MCMC
source("JAGS_R.R")

model.file <- "model.bug"

mod <- "
var
  n, x[n], y[n],
  y1,
  tau.x, sigma.x,
  tau.y, sigma.y;

model {
  x[1] ~ dnorm(y1, tau.x);
  y[1] ~ dnorm(x[1], tau.y);

  for (i in 2:length(y)) {
    x[i] ~ dnorm(x[i - 1], tau.x);
    y[i] ~ dnorm(x[i], tau.y);
  }

  tau.x ~ dgamma(1.0E-3, 1.0E-3);
  tau.y ~ dgamma(1.0E-3, 1.0E-3);
  sigma.x <- 1 / sqrt(tau.x);
  sigma.y <- 1 / sqrt(tau.y);
}
"

write(mod, file = model.file)

n.chains <- 3
data <- list(y = as.vector(Nile), n = length(Nile),
             y1 = Nile[1])
init1 <- list(.RNG.name = "base::Mersenne-Twister",
              .RNG.seed = 314,
              tau.x = 1, tau.y = 100)
init2 <- list(.RNG.name = "base::Mersenne-Twister",
              .RNG.seed = 3141,
              tau.x = 100, tau.y = 0.01)
init3 <- list(.RNG.name = "base::Mersenne-Twister",
              .RNG.seed = 31415,
              tau.x = 0.01, tau.y = 1)
inits <- list(init1, init2, init3)
var = c("x", "sigma.x", "sigma.y")

post <- jags.run(data = data, inits = inits,
                 parameters.to.save = var,
                 model.file = model.file,
                 n.chains = n.chains,
                 n.iter = 35000, n.burnin = 15000,
                 n.thin = 10)

gelman.diag(post)

summary(post[, (length(Nile) + 1):(length(Nile) + 2)])

y.mean <- apply(rbind(post[[1]][, 1:length(Nile)],
                      post[[2]][, 1:length(Nile)],
                      post[[3]][, 1:length(Nile)]), 2, mean)

グラフを重ねてみる。

## plot
s <- start(Nile)[1]
e <- end(Nile)[1]
plot(Nile, type ='o', las = 1)
lines(dropFirst(smoothNile$s), col = "red")
lines(s:e, y.mean, col = "blue")

タグ：MCMC R jags dlm