Как сохранить объект модели JAGS в R?

Я использую пакет rjags для выполнения MCMC в R и хочу сохранить вывод функции 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 для дальнейшего использования?

Примечание. Вопрос задавался ранее, но ОП был не очень конкретно о проблеме.


r jags
person sieste    schedule 21.08.2014    source источник

Ответы (1)


arrow_upward
7
arrow_downward

Может быть, я совершенно не в курсе того, что вы действительно хотите сделать, но я бы настроил такую ​​​​модель jags, используя R2jags вместо rjags (просто что-то вроде другой оболочки):

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 в список mcmc, который может обрабатываться кодом пакета.

Удачи!

P.S. Вот второй ответ:

После повторного просмотра вашего кода единственное, что после загрузки jags-объекта отсутствует для получения желаемого поведения, это:

jags$recompile()
coda.samples(model=jags,n.iter=1,variable.names="mu")

Но если вы действительно просто хотите использовать уже полученные апостериорные образцы или, может быть, просто хотите обновить цепочки для большего количества итераций, вы также можете использовать процедуру R2jags.

person Woosah    schedule 22.08.2014
comment
Ваш второй ответ отлично решил проблему. Большое спасибо. В своем первоначальном вопросе я не указал, что хочу восстановить сгоревшую цепь Маркова после запуска нового сеанса R. Я отредактировал вопрос соответственно. - person sieste; 22.08.2014