我正在使用该软件包rjags
在R中进行MCMC,我想保存该函数的输出以jags.model
供以后在另一个R会话中使用。
这是一个正态分布均值的简单示例:
library(rjags)
N <- 1000
x <- rnorm(N, 0, 5)
model.str <- 'model {for (i in 1:N) {
x[i] ~ dnorm(mu, 5)}
mu ~ dnorm(0, .0001)}'
jags <- jags.model(textConnection(model.str), data = list(x = x, N = N))
update(jags, 1000)
我可以生成这样的样本mu
:
coda.samples(model=jags,n.iter=1,variable.names="mu")
# [[1]]
# Markov Chain Monte Carlo (MCMC) output:
# Start = 2001
# End = 2001
# Thinning interval = 1
# mu
# [1,] 0.2312028
#
# attr(,"class")
# [1] "mcmc.list"
现在,我想保存模型对象jags
以供以后在新的R会话中使用,这样就不必再次初始化并在“马尔可夫链”中刻录:
save(file="/tmp/jags.Rdata", list="jags")
quit()
但是,在开始新的R会话并重新加载模型后,我收到一条错误消息,必须重新编译JAGS模型:
load("/tmp/jags.Rdata")
coda.samples(model=jags,n.iter=1,variable.names="mu")
# Error in model$iter() : JAGS model must be recompiled
这是为什么?如何将对象保存jags
在R中以供以后使用?
也许我对您真正想做的事情完全不了解,但是我会使用R2jags而不是rjags(就像其他包装器一样)来建立一个像这样的jags模型:
library(R2jags)
N <- 1000
x <- rnorm(N, 0, 5)
sink("test.txt")
cat("
model{
for (i in 1:N) {
x[i] ~ dnorm(mu, 5)
}
mu ~ dnorm(0, .0001)
}
",fill = TRUE)
sink()
inits <- function() {
list(
mu = dnorm(1, 0, 0.01))
}
params <- c("mu")
chains <- 3
iter <- 1000
jags1 <- jags(model.file = "test.txt", data = list(x = x, N = N),
parameters.to.save = params, inits = inits,
n.chains = chains, n.iter = iter, n.burnin=floor(iter/2),
n.thin = ifelse(floor(iter/100) < 1, 1, floor(iter/100)))
jags2 <- update(jags1, 10000)
jags2
plot(jags2)
traceplot(jags2)
jags2.mcmc <- as.mcmc(jags2)
结果没有差异,我喜欢这个过程,因为它更多地是我使用winbugs的方式,所以...
代码的最后一行将jags2-object转换为mcmc-list,可以通过包coda对其进行处理。
祝你好运!
PS这是第二个答案:
再次查看代码后,加载无法获得所需行为的jags对象后,唯一的事情是:
jags$recompile()
coda.samples(model=jags,n.iter=1,variable.names="mu")
但是,如果您真的只想使用已经获得的后验样本,或者只是想更新链以进行更多迭代,则也可以使用R2jags-procedure。
本文收集自互联网,转载请注明来源。
如有侵权,请联系 [email protected] 删除。
我来说两句