So-net無料ブログ作成

1次元CARと状態空間モデルとの等価性 [統計]

前回の記事についてbaibai先生から、dlmと尤度がガウスの1次元CARで結果を比較してみるとよいかもといわれましたので、やってみました。

前回のJAGSのモデルのうち、観測モデルの部分を正規分布にしたものと、dlmの結果とを比較してみます。dlmのシステムモデルはローカルレベルモデルです。

Rコード

Y <- c( 0,  3,  2,  5,  6, 16,  8, 14, 11, 10,
       17, 19, 14, 19, 19, 18, 15, 13, 13,  9,
       11, 15, 18, 12, 11, 17, 14, 16, 15,  9,
        6, 15, 10, 11, 14,  7, 14, 14, 13, 17,
        8,  7, 10,  4,  5,  5,  7,  4,  3,  1)
N.site <- length(Y)

## dlm
library(dlm)
build1 <- function(theta) {
  dlmModPoly(order = 1, dV = theta[1], dW = theta[2])
}
dlm.fit1 <- dlmMLE(Y, build = build1, parm = c(1, 1), lower = rep(1e-8, 2))
Y.dlm <- dlmSmooth(Y, build1(dlm.fit1$par))$s[-1]

## JAGS: Gauss
library(rjags)
model <- jags.model("jags_car2_bug.txt",
                    data = list(Y = Y, N.site = N.site),
                    n.chains = 3, n.adapt = 5000)
fit <- coda.samples(model,
                    variable.names = c("r", "sigma"),
                    n.iter = 10000, thin = 10)
Y.jags <- summary(fit[, grep("r", varnames(fit))])$statistics[, "Mean"]

library(ggplot2)
df <- data.frame(site = seq_along(Y),
                 Y = Y)
df2 <- data.frame(site = rep(seq_along(Y), 2),
                  method = rep(c("dlm", "JAGS"), each = N.site),
                  Y = c(Y.dlm,Y.hat))
ggplot(df) +
  geom_point(aes(x = site, y = Y), size = 3) +
  geom_line(data = df2, aes(x = site, y = Y, colour = method), alpha = 0.5, size = 2)

JAGSコード

var
  N.site,
  Y[N.site],
  r[N.site],
  sigma[2],
  tau[2];
model {
  for (i in 1:N.site) {
    Y[i] ~ dnorm(r[i], tau[1]);
  }
  for (i in 2:N.site) {
    r[i] ~ dnorm(r[i - 1], tau[2]);
  }
  r[1] ~ dnorm(0, 1.0E-4);
  for (i in 1:2) {
    tau[i] <- 1 / (sigma[i] * sigma[i]);
    sigma[i] ~ dunif(0, 1.0E+4);
  }
}

結果

Rplot002.png

dlmでは最尤推定値を使用し、JAGSでは事後平均を使用したのですが、だいたいぴったりあっています。

では、OpenBUGSでcar.normalを使用した場合と比較するとどうか、ということでやってみました。JAGSは今回は観測モデルをポアソン分布とし(前回のモデルと同じなのでコード省略)、ポアソン分布の状態空間モデルを使えるKFASもいれてみました。

Rコード(OpenBUGSの部分は、自分の環境にあわせてWine経由で使用しているので、適宜修正してください)

Y <- c( 0,  3,  2,  5,  6, 16,  8, 14, 11, 10,
       17, 19, 14, 19, 19, 18, 15, 13, 13,  9,
       11, 15, 18, 12, 11, 17, 14, 16, 15,  9,
        6, 15, 10, 11, 14,  7, 14, 14, 13, 17,
        8,  7, 10,  4,  5,  5,  7,  4,  3,  1)
N.site <- length(Y)

## JAGS
library(rjags)
model <- jags.model("jags_car_bug.txt",
                    data = list(Y = Y, N.site = N.site),
                    n.chains = 3, n.adapt = 5000)
fit <- coda.samples(model,
                    variable.names = c("r", "sigma"),
                    n.iter = 10000, thin = 10)

Y.jags <- exp(summary(fit[, grep("r", varnames(fit))])$statistics[, "Mean"])

## OpenBUGS
library(R2OpenBUGS)
library(coda)

WINE <- "/Applications/Wine.app/Contents/Resources/bin/wine"
WINEPATH <- "/Applications/Wine.app/Contents/Resources/bin/winepath"
OpenBUGS <- paste(Sys.getenv("HOME"),
                  ".wine/drive_c/Program Files/OpenBUGS/OpenBUGS323/OpenBUGS.exe",
                  sep = "/")
adj <- c(2, c(sapply(2:49, function(i) c(i - 1, i + 1))), 49)
weights <- rep(1, length(adj))
num <- c(1, rep(2, 48), 1)
inits <- list(list(S = rep(0, 50), beta = 1, sigma = 1),
              list(S = rep(1, 50), beta = 2, sigma = 2),
              list(S = rep(2, 50), beta = 3, sigma = 3))

bugs.fit <- bugs(data = list(N.site = length(Y), Y = Y,
                             Adj = adj, Weights = weights, Num = num),
                 inits = inits,
                 parameters.to.save = c("beta", "r", "sigma"),
                 model.file = "CAR_bug.txt",
                 n.chains = 3, n.iter = 11000, n.burnin = 1000, n.thin = 10,
                 OpenBUGS.pgm = OpenBUGS,
                 WINE = WINE, WINEPATH = WINEPATH,
                 useWINE = TRUE)
Y.bugs <- exp(bugs.fit$mean$beta + bugs.fit$mean$r)

## dlm
library(dlm)

## 1階差分
build1 <- function(theta) {
  dlmModPoly(order = 1, dV = theta[1], dW = theta[2])
}
dlm.fit1 <- dlmMLE(Y, build = build1, parm = c(1, 1), lower = rep(1e-8, 2))
Y.dlm <- dlmSmooth(Y, build1(dlm.fit1$par))$s[-1]

## KFAS
library(KFAS)

## 1階差分
kfas.model1 <- SSModel(Y ~ SSMtrend(degree = 1, Q = list(matrix(NA))),
                       distribution = "poisson")
kfas.fit1 <- fitSSM(kfas.model1, inits = c(1, 1))
kfas.out1 <- KFS(kfas.fit1$model, smoothing = c("state"))
kfas.sig1 <- signal(kfas.out1, states = c("all"))

Y.kfas <- exp(kfas.sig1$signal)

## ggplot
library(ggplot2)
df <- data.frame(site = seq_along(Y),
                 Y = Y)
df2 <- data.frame(site = rep(seq_along(Y), 4),
                  method = rep(c("JAGS", "OpenBUGS", "dlm", "KFAS"),
                               each = length(Y)),
                  Y.hat = c(Y.jags, Y.bugs, Y.dlm, Y.kfas))
ggplot(df) +
  geom_point(aes(x = site, y = Y), size = 3) +
  geom_line(data = df2, aes(x = site, y = Y.hat, colour = method),
            size = 2, alpha = 0.5)

OpenBUGSのコード

model {
  for (i in 1:N.site) {
    Y[i] ~ dpois(lambda[i]);
    log(lambda[i]) <- beta + r[i];
  }
  r[1:N.site] ~ car.normal(Adj[], Weights[], Num[], tau);
  beta ~ dnorm(0, 1.0E-4);
  tau <- 1 / (sigma * sigma);
  sigma ~ dunif(0, 1.0E+4);
}

結果です。

Rplot001.png

OpenBUGSでcar.normalを使用した結果と、JAGSおよびKFASで観測モデルがポアソンの状態空間モデルを使用した結果とはほとんど重なっています。さらに、dlmの結果もわりと重なっていました。

ということで、1次元CARと状態空間モデルとの対応がわかる結果になりました。


タグ:R jags OpenBUGS
nice!(3)  コメント(0)  トラックバック(0) 
共通テーマ:日記・雑感

nice! 3

コメント 0

コメントを書く

お名前:
URL:
コメント:
画像認証:
下の画像に表示されている文字を入力してください。

Facebook コメント

トラックバック 0