Управление объектом 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 секунды.