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

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

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

Однако я не смог понять, как получить апостериорные рисунки из 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"

Привет, Джейкоб


r mcmc data-manipulation jags
person Jacob Socolar    schedule 15.11.2015    source источник
comment
В rjags может быть метод, но я, кажется, помню, что вы можете сделать do.call(rbind.data.frame, foldD_samples). Возможно, жирнее эффективнее использовать data.table::rbindlist   -  person user20650    schedule 15.11.2015
comment
PS вы можете применить code.samples к списку, не принуждая к кадру данных   -  person user20650    schedule 15.11.2015
comment
Спасибо пользователю 20650! do.call(rbind.data.frame, foldD_samples) работал хорошо. Я был бы рад принять это как ответ, если бы он был опубликован как таковой. data.table::rbindlist не принимает объекты mcmc.list в качестве входных данных. Также обратите внимание на предполагаемый код опечатки для коды в постскриптуме.   -  person Jacob Socolar    schedule 15.11.2015
comment
rbindlist(lapply(foldD_samples,as.data.frame)) может сработать... если вам важна эффективность.   -  person Ben Bolker    schedule 15.11.2015

Ответы (2)


arrow_upward
6
arrow_downward

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

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

or

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)
person user20650    schedule 15.11.2015

arrow_upward
2
arrow_downward

Ответ, данный пользователем 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 цепочек, всего 50 000 строк, метод 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 секунды за прошедшее время.

person prince_of_pears    schedule 16.07.2017