Управление объектом mcmc.list в R

Я использовал JAGS, вызванный через rjags, для создания объекта mcmc.list foldD_samples, который содержит мониторы трассировки для большого количества стохастических узлов (>800 узлов).

Теперь я хотел бы использовать R для вычисления довольно сложной, скалярной функции этих узлов и записать вывод в объект mcmc, чтобы я мог использовать coda для суммирования апостериорной и прогностической диагностики сходимости.

Тем не менее, я не смог выяснить, как получить задний рисунок из foldD_samples в кадр данных. Любая помощь высоко ценится.

Вот структура списка mcmc.list:

str(foldD_samples)
List of 3
 $ : mcmc [1:5000, 1:821] -0.667 -0.197 -0.302 -0.204 -0.394 ...
  ..- attr(*, "dimnames")=List of 2
  .. ..$ : NULL
  .. ..$ : chr [1:821] "beta0" "beta1" "beta2" "dtau" ...
  ..- attr(*, "mcpar")= num [1:3] 4100 504000 100
 $ : mcmc [1:5000, 1:821] -0.686 -0.385 -0.53 -0.457 -0.519 ...
  ..- attr(*, "dimnames")=List of 2
  .. ..$ : NULL
  .. ..$ : chr [1:821] "beta0" "beta1" "beta2" "dtau" ...
  ..- attr(*, "mcpar")= num [1:3] 4100 504000 100
 $ : mcmc [1:5000, 1:821] -0.492 -0.679 -0.299 -0.429 -0.421 ...
  ..- attr(*, "dimnames")=List of 2
  .. ..$ : NULL
  .. ..$ : chr [1:821] "beta0" "beta1" "beta2" "dtau" ...
  ..- attr(*, "mcpar")= num [1:3] 4100 504000 100
 - attr(*, "class")= chr "mcmc.list"

Ура, Джейкоб

2 ответа

Решение

Как это list Структуру вы можете использовать любой из этих методов, чтобы связать матрицы вместе.

do.call(rbind.data.frame, foldD_samples)

или же

rbindlist(lapply(foldD_samples, as.data.frame)) # thanks to BenBolker

Мве

library(rjags)
library(coda)
library(data.table)

mod <- textConnection("model {
  A ~ dnorm(0, 1)
  B ~ dnorm(0, 1)
}")

# evaluate
mod <- jags.model(mod, n.chains = 4, n.adapt = 50000) 
pos <- coda.samples(mod,  c("A", "B"),  10000)

out <- do.call(rbind.data.frame, pos)
out2 <- rbindlist(lapply(pos, as.data.frame))
all.equal(out, out2, check.attributes=FALSE)

Ответ пользователя 20650, безусловно, будет работать, но он может быть довольно медленным. Также обратите внимание, что на момент написания этой статьи rbind_list() устарел вместо bind_rows().

Что-то, что я написал для своих собственных целей, преобразует mcmc.list в "длинный формат" data.frame. На моей машине это примерно в 4-7 раз быстрее, чем описанные выше методы, и добавляет два дополнительных столбца: один для номера цепи и один для номера шага.

parameter_names <- varnames(mcmc_list)
saved_steps <- as.integer(row.names(mcmc_list[[1]]))
out <- data.frame("chain" = factor(rep(1 : length(mcmc_list), each = length(saved_steps))),
                  "step" = rep(saved_steps, length(mcmc_list)) )
for (param in parameter_names) {
    out[param] <- NA
}
for (a_chain in 1 : length(mcmc_list)) {
    out[out$chain == a_chain, parameter_names ] <- as.data.frame(mcmc_list[[a_chain]])
}
return(out)

Работа с объектом mcmc.list из 3 цепочек, всего 50000 строк, метод rbind_list/do.call: 0,71 с, истекшее время мой метод: 0,12 с, истекшее время

Изменить: некоторое дальнейшее чтение кода в библиотеке "coda" показывает, что as.matrix() намного быстрее.

parameter_names <- varnames(mcmc_list)
saved_steps <- as.integer(row.names(mcmc_list[[1]]))
out <- data.frame("chain" = factor(rep(1 : length(mcmc_list), each = length(saved_steps))),
                  "step" = rep(saved_steps, length(mcmc_list)) )
out <- cbind(out, as.data.frame(as.matrix(chain_samples)))

Занимает более 0,03 секунды.

Другие вопросы по тегам