R: mcmc.listをlatticeで表示

R: mcmc.listをggplot2で表示につづいて、こんどはlatticeグラフィックスで表示してみる。買ったばかりの本を参考に。



コード

library(lattice)
library(latticeExtra)
library(rjags)

## data
K <- 10
X <- c(0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10)
Y <- c(1, 2, 2, 6, 4, 5, 8, 9, 9, 9, 10)
N <- length(X)

model.text <- "
model {
  for (i in 1:N) {
    Y[i] ~ dbin(p[i], K);
    logit(p[i]) <- beta + beta.x * X[i];
  }
  
  # priors
  beta ~ dnorm(0.0, 1.0E-6);
  beta.x ~ dnorm(0.0, 1.0E-6);
}
"
cat(model.text, file = "model.bug")

model <- jags.model("model.bug",
                    data = list(K = K, X = X, Y = Y, N = N),
                    n.chains =3,
                    n.adapt = 1000)
samp <- coda.samples(model,
                     variable.names = c("beta", "beta.x"),
                     n.iter = 20000,
                     thin = 10)

x <- rep(1:niter(samp))
beta <- beta.x <- vector("list", nchain(samp))
for (i in 1:nchain(samp)) {
  beta[[i]] <- samp[[i]][, "beta"]
  beta.x[[i]] <- samp[[i]][, "beta.x"]
}

pdf("mcmc_lattice.pdf", width = 256/72, height = 256/72,
    family = "Helvetica", pointsize = 9)
## trace
xyplot(beta[[1]] + beta[[2]] + beta[[3]] ~ x, type = "l",
       xlab = "iterations", ylab = "beta")
xyplot(beta.x[[1]] + beta.x[[2]] + beta.x[[3]] ~ x, type = "l",
       xlab = "iterations", ylab = "beta.x")

## density
densityplot(~c(beta[[1]], beta[[2]], beta[[3]]),
            plot.points = "rug",
            xlab = "beta", ylab = "density")
densityplot(~c(beta.x[[1]], beta.x[[2]], beta.x[[3]]),
            plot.points = "rug",
            xlab = "beta.x", ylab = "density")
dev.off()

タグ：R Lattice