読者です 読者をやめる 読者になる 読者になる

StatModeling Memorandum

StanとRとPythonでベイズ統計モデリングします. たまに書評.

JAGSで2次元のマルコフ場モデル

WinBUGS, JAGS, Stanの三姉妹を等しく愛する僕としては何使っても同時分布を自由に使えるんだよ、ということで2次元データのマルコフ場モデルをやりました。

まず1次元データに対するCAR modelは、Itoさんのこの記事この記事を参照。Rパッケージの{dlm}{KFAS}を使った場合との比較があり大変参考になります。

本題の2次元の場合。データは前の記事の後半のものを使います。

JAGSで自由に尤度を設定したい場合はzero trickを使います(例えばこの記事参照)。

それではJAGSコード行きます。まずは1次(1階差分)のマルコフ場モデルから。

data {
   for (i in 1:N){
      for (j in 2:N){
         z.h[i,j-1] <- 0
      }
   }
   for (j in 1:N){
      for (i in 2:N){
         z.v[i-1,j] <- 0
      }
   }
}

model {
   for (i in 1:N){
      for (j in 2:N){
         z.h[i,j-1] ~ dpois(neg.log.lik.h[i,j-1])
         neg.log.lik.h[i,j-1] <- log(s.r) + 0.5*pow((r[i,j] - r[i,j-1])/s.r, 2) + Const
      }
   }
   for (j in 1:N){
      for (i in 2:N){
         z.v[i-1,j] ~ dpois(neg.log.lik.v[i-1,j])
         neg.log.lik.v[i-1,j] <- log(s.r) + 0.5*pow((r[i,j] - r[i-1,j])/s.r, 2) + Const
      }
   }
   s.r ~ dunif(0, 1000)

   for (i in 1:N){
      for (j in 1:N){
         r[i,j] ~ dnorm(0, 1.0E-4)
         Z[i,j] ~ dnorm(r[i,j], tau.z)
      }
   }
   tau.z <- pow(s.z, -2)
   s.z ~ dunif(0, 1000)
}
  • 31行目: rはパラメータということを明示するために平らな事前分布でよいのでこうやって書いておく必要があります。

2次(2階差分)のマルコフ場モデルのJAGSコードは以下のようになります。

data {
   for (i in 1:N){
      for (j in 3:N){
         z.h[i,j-2] <- 0
      }
   }
   for (j in 1:N){
      for (i in 3:N){
         z.v[i-2,j] <- 0
      }
   }
}

model {
   for (i in 1:N){
      for (j in 3:N){
         z.h[i,j-2] ~ dpois(neg.log.lik.h[i,j-2])
         neg.log.lik.h[i,j-2] <- log(s.r) + 0.5*pow((r[i,j] - 2*r[i,j-1] + r[i,j-2])/s.r, 2) + Const
      }
   }
   for (j in 1:N){
      for (i in 3:N){
         z.v[i-2,j] ~ dpois(neg.log.lik.v[i-2,j])
         neg.log.lik.v[i-2,j] <- log(s.r) + 0.5*pow((r[i,j] - 2*r[i-1,j] + r[i-2,j])/s.r, 2) + Const
      }
   }
   s.r ~ dunif(0, 1000)

   for (i in 1:N){
      for (j in 1:N){
         r[i,j] ~ dnorm(0, 1.0E-4)
         Z[i,j] ~ dnorm(r[i,j], tau.z)
      }
   }
   tau.z <- pow(s.z, -2)
   s.z ~ dunif(0, 1000)
}
  • 18,24行目: 特に変化があるのはここだけです。

1次のマルコフ場モデルの場合の、キックするRコードは以下になります。

library(rjags)

N <- 30
x <- seq(-10, 10, length=N)
y <- x
f <- function(x, y) { r <- sqrt(x^2+y^2); 10*sin(r)/r }
z <- outer(x, y, f)
set.seed(123)
noise <- matrix(rnorm(N*N, 0, 0.5), nrow=N, ncol=N)
Z <- z + noise

data <- list(N=N, Z=Z, Const=100)
inits <- list(r=Z, s.r=1, s.z=1)

m <- jags.model('model/model1.bugs', data, inits, n.chains=3, n.adapt=1000)
update(m, 3000)
post.list <- coda.samples(m, c('r', 's.r', 's.z'), n.iter=10000, thin=10)

library(R2WinBUGS)
mcmc.list2bugs <- function(mcmc.list){
   b1 <- mcmc.list[[1]]
   m1 <- as.matrix(b1)
   mall <- matrix(numeric(0), 0, ncol(m1))
   n.chains <- length(mcmc.list)
   for (i in 1:n.chains) {
      mall <- rbind(mall, as.matrix(mcmc.list[[i]]))
   }
   sims.array <- array(mall, dim = c(nrow(m1), n.chains, ncol(m1)))
   dimnames(sims.array) <- list(NULL, NULL, colnames(m1))
   mcpar <- attr(b1, "mcpar")
   as.bugs.array(
      sims.array = sims.array,
      model.file = NULL,
      program = NULL,
      DIC = FALSE,
      DICOutput = NULL,
      n.iter = mcpar[2],
      n.burnin = mcpar[1] - mcpar[3],
      n.thin = mcpar[3]
   )
}
post.bugs <- mcmc.list2bugs(post.list)
post.mcmc <- as.mcmc(post.bugs$sims.matrix)
save.image('result.RData')


r_median <- t(sapply(1:N, function(i){
   tmp <- post.mcmc[ , sprintf('r[%d,%d]', i, 1:N)]
   apply(tmp, 2, median)
}))
png('persp-M1-r-median.png', pointsize=24, width=800, height=600)
persp(x, y, r_median, theta=30, phi=30, expand=0.5, col=rainbow(N), border=NA)
dev.off()
  • 13行目:前回同様rの初期値をかなりちゃんと設定しないとうまく収束しませんでした。 なお計算時間はSurface Pro 3(core i5)で2次のマルコフ場モデルが3chainで5分ぐらいでした。意外と早い。

結果はパラメータrのMCMCサンプルの中央値だけとってきて、3次元でplotしたものが以下になります。左図が1次・右図が2次。Stanの時と結果がほぼ同じであることが分かります。enjoy!

20150223_persp-M1-r-median.png20150223_persp-M2-r-median.png