JAGS: nested indexing vs matrix

混合モデルを階層ベイズモデルで解くときに、データを行列にするか、nested indexingを使うか、という問題。どちらが速いか、JAGSで試してみた。

こういうコードで。test1.Rはnested indexing。

test1.R

## model 1

library(rjags)

set.seed(11)

i.logit <- function(x) exp(x)/(1+exp(x))

b <- 1:10
m <- length(b)
e <- rnorm(m, 0, 0.5)
n <- 21
x <- seq(-1, 1, length = n)
data <- expand.grid(x = x, b = b)
data$y <- rbinom(m * n, 1,
                 i.logit(0.8 + 1.5 * data$x + e[data$b]))

model1 <- jags.model("model1.bug",
                     data = list(K = n * m,
                                 M = m,
                                 y = data$y,
                                 x = data$x,
                                 b = data$b),
                      nchain = 3, n.adapt = 5000)
post1 <- coda.samples(model1,
                      variable = c("beta", "beta.x", "sigma"),
                      n.iter = 10000, thin = 10)
summary(post1)

model1.bug

model {
  for (i in 1:K) {
    y[i] ~ dbern(p[i]);
    logit(p[i]) <- beta + beta.x * x[i] + e[b[i]];
  }
  for (j in 1:M) {
    e[j] ~ dnorm(0, tau);
  }

  beta ~ dnorm(0.0, 1.0E-4);
  beta.x ~ dnorm(0.0, 1.0E-4);
  tau ~ dgamma(1.0E-2, 1.0E-2);
  sigma <- 1/sqrt(tau);
}

続いて行列で。

test2.R

## model 2

library(rjags)

set.seed(11)

i.logit <- function(x) exp(x)/(1+exp(x))

b <- 1:10
m <- length(b)
e <- rnorm(m, 0, 0.5)
n <- 21
x <- seq(-1, 1, length = n)
data <- expand.grid(x = x, b = b)
data$y <- rbinom(m * n, 1,
                 i.logit(0.8 + 1.5 * data$x + e[data$b]))

model2 <- jags.model("model2.bug",
                     data = list(N = n,
                                 M = m,
                                 y = matrix(data$y, ncol = m),
                                 x = data$x),
                      nchain = 3, n.adapt = 5000)
post2 <- coda.samples(model2,
                      variable = c("beta", "beta.x", "sigma"),
                      n.iter = 10000, thin = 10)
summary(post2)

model2.bug

model {
  for (j in 1:M) {
    for (i in 1:N) {
      y[i,j] ~ dbern(p[i,j]);
      logit(p[i,j]) <- beta + beta.x * x[i] + e[j];
    }
    e[j] ~ dnorm(0, tau);
  }

  beta ~ dnorm(0.0, 1.0E-4);
  beta.x ~ dnorm(0.0, 1.0E-4);
  tau ~ dgamma(1.0E-2, 1.0E-2);
  sigma <- 1/sqrt(tau);
}

で、結果はというと

test1

Iterations = 5010:15000
Thinning interval = 10 
Number of chains = 3 
Sample size per chain = 1000 

1. Empirical mean and standard deviation for each variable,
   plus standard error of the mean:

         Mean     SD Naive SE Time-series SE
beta   0.8165 0.2120 0.003871       0.004399
beta.x 1.9319 0.3171 0.005789       0.005553
sigma  0.3291 0.2100 0.003835       0.005376

2. Quantiles for each variable:

          2.5%    25%    50%    75% 97.5%
beta   0.41311 0.6768 0.8106 0.9485 1.242
beta.x 1.33061 1.7119 1.9208 2.1378 2.588
sigma  0.08182 0.1754 0.2809 0.4279 0.861

> 
> proc.time()
   ユーザ   システム       経過  
    49.206      0.195     50.195 

test2

Iterations = 5010:15000
Thinning interval = 10 
Number of chains = 3 
Sample size per chain = 1000 

1. Empirical mean and standard deviation for each variable,
   plus standard error of the mean:

         Mean     SD Naive SE Time-series SE
beta   0.8171 0.2181 0.003981       0.004316
beta.x 1.9369 0.3137 0.005728       0.006621
sigma  0.3201 0.2197 0.004011       0.005987

2. Quantiles for each variable:

          2.5%    25%    50%    75%  97.5%
beta   0.39813 0.6766 0.8119 0.9499 1.2579
beta.x 1.34165 1.7256 1.9271 2.1345 2.5853
sigma  0.07614 0.1657 0.2607 0.4177 0.8955

> 
> proc.time()
   ユーザ   システム       経過  
    49.212      0.294     50.033 

どちらでもあまり変わらないみたい。

タグ：jags