我想实现以下模型JAGS为WinBUGS软件写:翻译WinBUGS软件模型JAGS(用R)
model {
for (i in 1:N) {
wtp[i] ~ dweib(r[G[i]], mu[i])I(lower[i], upper[i])
mu[i] <- exp(beta[G[i]])
G[i] ~ dcat(P[])
}
P[1] ~ dunif(0.01, 0.99)
P[2] <- 1 - P[1]
r[1] ~ dunif(1, 10)
r[2] ~ dunif(0.1, 10)
beta[1] ~ dunif(0, 1000)
beta[2] ~ dunif(-1000, 0)
weibmed[1] <- pow(log(2) * exp(-beta[1]), 1/r[1])
weibmed[2] <- pow(log(2) * exp(-beta[2]), 1/r[2])
weibmed[3] <- pow(log(1/(1 - 0.5 + P[1])) * exp(-beta[2]), 1/r[2])
weibmean[1] <- pow(exp(-beta[1]), 1/r[1]) * exp(loggam((1 + r[1])/r[1]))
weibmean[2] <- pow(exp(-beta[2]), 1/r[2]) * exp(loggam((1 + r[2])/r[2]))
weibmean[3] <- P[1] * weibmean[1] + P[2] * weibmean[2]
}
我认为这将是简单的把它在JAGS要和:
library(rjags)
txt <- 'model {
for (i in 1:N) {
wtp[i] ~ dweib(r[G[i]], mu[i])T(lower[i], upper[i])
mu[i] <- exp(beta[G[i]])
G[i] ~ dcat(P[])
}
P[1] ~ dunif(0.01, 0.99)
P[2] <- 1 - P[1]
r[1] ~ dunif(1, 10)
r[2] ~ dunif(0.1, 10)
beta[1] ~ dunif(0, 1000)
beta[2] ~ dunif(-1000, 0)
weibmed[1] <- pow(log(2) * exp(-beta[1]), 1/r[1])
weibmed[2] <- pow(log(2) * exp(-beta[2]), 1/r[2])
weibmed[3] <- pow(log(1/(1 - 0.5 + P[1])) * exp(-beta[2]), 1/r[2])
weibmean[1] <- pow(exp(-beta[1]), 1/r[1]) * exp(loggam((1 + r[1])/r[1]))
weibmean[2] <- pow(exp(-beta[2]), 1/r[2]) * exp(loggam((1 + r[2])/r[2]))
weibmean[3] <- P[1] * weibmean[1] + P[2] * weibmean[2]
}'
set.seed(3.14159)
dat <- list(N = 1000, lower = rep(0, 1000), upper = runif(1000, 5, 200000))
ini <- list(P = c(0.4, NA), r = c(8.2, 1.2), beta = c(3.8, -6.5))
mod <- jags.model(
file = textConnection(txt),
data = dat,
inits = c(ini, .RNG.name = 'base::Mersenne-Twister', .RNG.seed = 314159),
n.chains = 1,
n.adapt = 100
)
sam.jags <- coda.samples(
model = mod,
variable.names = c('P', 'r', 'beta', 'weibmed', 'weibmean'),
n.iter = 400,
n.thin = 1
)
只需将I()
替换为T()
即可。这将产生coda.samples()
错误:
Error: Error in node weibmed[3]
Invalid parent values
如果我忽略的weibmed
监测和weibmean
然后coda.samples()
作品,但参数估计值:
Mean SD Naive SE Time-series SE
P[1] 0.4840704 0.2769491 0.01384746 0.01384746
P[2] 0.5159296 0.2769491 0.01384746 0.01384746
beta[1] 509.3614647 295.0860473 14.75430237 14.75430237
beta[2] -487.5362940 285.4126899 14.27063449 14.27063449
r[1] 5.2054730 2.6330434 0.13165217 0.13165217
r[2] 5.0478143 2.9480476 0.14740238 0.14740238
没有可比性那些我使用WinBUGS软件时得到:
library(R2WinBUGS)
sam.bugs <- bugs(
model.file = 'model.bug',
data = dat,
inits = list(ini),
parameters.to.save = c('P', 'r', 'beta'), #, 'weibmed', 'weibmean'),
n.chains = 1,
n.burnin = 100,
n.iter = 500,
n.thin = 1,
debug = F,
DIC = F,
bugs.seed = 314159
)
Inference for Bugs model at "3mixout2.bug", fit using WinBUGS,
1 chains, each with 500 iterations (first 100 discarded)
n.sims = 400 iterations saved
mean sd 2.5% 25% 50% 75% 97.5%
P[1] 0.4 0.0 0.3 0.4 0.4 0.4 0.4
P[2] 0.6 0.0 0.6 0.6 0.6 0.6 0.7
r[1] 7.2 0.8 6.2 6.6 6.9 7.7 9.3
r[2] 1.5 0.0 1.4 1.4 1.5 1.5 1.6
beta[1] 5.5 0.5 4.6 5.0 5.4 5.8 6.6
beta[2] -7.2 0.2 -7.5 -7.3 -7.2 -7.1 -6.9
有什么想法或建议吗?
感谢您指出了这一点。除了受监控的变量列表中的“weibmed”和“weibmean”,你是否还有理由相信WinBUGS和JAGS会为'P','r'和'beta'产生不同的估计值? –
@ feats-by-jake;对不起,我把旧电脑带走了BUGS环境。你能告诉我结果吗? (以及它是否通过'parameters.to.save'除了'weibmed'和'weibmean'之外改变了吗?) – cuttlefish44
在这两种情况下,我都同时排除了'weibmed'和'webmean',并在'coda.samples() '和'bugs()'。如果'weibmed'和'webmean'被添加回来,WinBUGS输出不会改变。 –