trapの呪いから抜けられない

どうにか、zipモデルをwinbugsコードで書くことはできたのですが、trapから抜けられそうな気配が全くありません。
ちなみにcongdon先生バージョンzipモデルはこういう風に書いてみた。
model {
# Likelihood
# Loop for Plot
for (i in 1:n) {
log(output[i]) <- alpha + beta*elv[i] + w[i]
deer[i] ~ dpois(mu.i[i])
mu.i[i] <- m[i,index[i]]
m[i,1] <- 0
m[i,2] <- output[i]
index[i] ~ dcat(theta[i,1:2])
theta[i,1] <- p[i]
theta[i,2] <- 1 – p[i]
logit(p[i]) <- betaz *elv[i]
mu[i] <- 0
}
w[1:n] ~ spatial.exp(mu[],x[],y[],tau,phi,1)
# Parameters and hyper parameters
alpha ~ dnorm(0,0.1)
beta ~ dnorm(0,0.1)
betaz ~ dnorm(0,0.1)
tau ~ dgamma(0.001, 0.1)
phi ~ dunif(0.05, 20)
sigma <- 1/tau
}
kuboさんバージョンzipモデルはこーゆー感じ。
model {
# Likelihood
# Loop for Plot
for (i in 1:n) {
deer[i] ~ dpois(out[i])
out[i] <- (1 – zero[i]) *p.n[i]
log(p.n[i]) <- alpha
zero[i] ~ dbern(p.z[i])
#logit(p.z[i]) <- betaz * elv[i] + w[i]
logit(p.z[i]) <- betaz + w[i]
mu[i] <- 0
}
w[1:n] ~ spatial.exp(mu[],x[],y[],tau,phi,1)
# Parameters and hyper parameters
alpha ~ dnorm(0,0.1)
beta ~ dnorm(0,0.1)
betaz ~ dnorm(0,0.1)
tau ~ dgamma(0.001, 0.1)
phi ~ dunif(0.05, 20)
sigma <- 1/tau
}
まずは、変数をscale()したら怒られたので、ベイジアンクリギングポアソンモデルでやったときのパラメータ値を初期値にぶっ込んでもみたが駄目。
次には暗黒面の奥義、事前分布操作にまで手を出したもののこれも駄目。
次に、変数減らして定数項のみ、つまりモデル自体何の意味もなく、ただMCMC回したいだけ、とゆー本末転倒なコードを流すと20分ほど回ってまたしてもtrap・・・・
あ゛~~~ぐあ~~~っとPCぶっ壊したい衝動に駆られながら、気持ちを落ち着けて冷静に考えてみても次なる手が浮かばない。
まぢ困った。
-
前の記事
花火見てきた 2008.08.08
-
次の記事
marryとかに関する何かのこもごも 2008.08.13