BUGS coding ばっどのうはう
……
というか，
いまいち見とおしのよくない書きかたなんだけど，
他にやりようがなさそうな一例．
X と Y というデータ (vector) があって，
X の中には NA (欠測) がまじっている．
X が NA である場合とそうでない場合で Y の分布を変えるには，
どうしたらよいか?
いろいろやりかたがあるんだけど，
一番素朴には下のように coding できる．
source("http://hosho.ees.hokudai.ac.jp/~kubo/ce/r/R2WBwrapper.R")
clear.data.param()
n <- 100
set.data("Y", rpois(n, lambda = 3.5)) # 架空データ
X <- sample(c(1:3, NA), n, replace = TRUE) # NA まじりの架空データ
X.isNA <- is.na(X)
set.data("V.Na", c(1:n)[X.isNA]) # NA であるデータ番号
set.data("V.NotNa", c(1:n)[!X.isNA]) # NA でないデータ番号
set.data("Length.V.Na", length(V.Na))
set.data("Length.V.NotNa", length(V.NotNa))
modelf <- function() # ここから R ソースファイル内に BUGS code をうめこむ
{
# X == NA data
for (i in 1:Length.V.Na) {
Y[V.Na[i]] ~ dpois(mean1)
}
# X != NA data
for (i in 1:Length.V.NotNa) {
Y[V.NotNa[i]] ~ dpois(mean2)
}
# parameters
log(mean1) <- log.mean1
log.mean1 ~ dnorm(0.0, 1.0E-4)
log(mean2) <- log.mean2
log.mean2 ~ dnorm(0.0, 1.0E-4)
}
write.model(modelf, "model.bug.txt") # ここから下はまた R
set.param("log.mean1", 0)
set.param("log.mean2", 0)
post.bugs <- call.bugs(
file = "model.bug.txt",
n.chains = 3,
debug = FALSE,
n.iter = 3000, n.burnin = 500, n.thin = 5
)
plot(post.bugs)
