Есть ли способ сохранить скорость принятия алгоритма MCMC, используемого MCMClogit?

Пакет R MCMCpack предлагает байесовскую логистическую регрессию через MCMClogit. Эта функция выводит скорость принятия алгоритма MCMC (Метрополис-Гастингс), если verbose=TRUE, но, похоже, не возвращает его в объекте mcmc. Есть ли способ получить доступ к этой информации в объекте?

Для проверки используйте:

    library(MCMCpack)
    set.seed(12345) 
    n = 1000 
    x = rnorm(n) 
    y = rbinom(n,1,1/(1+exp(-(1 + x))))
    m = MCMClogit(y ~ x, burnin = 5000, mcmc = 1000,
                         tune = 1.3, B0 = 0, verbose = TRUE)

Который печатает коэффициент принятия 0.45533, но я не нахожу эту информацию в names(m), возвращающем NULL или names(attributes(m)), который возвращает

[1] "dim"      "mcpar"    "class"    "dimnames" "title"    "y"        "call"

В файле справки указано, что пакет coda позволяет суммировать информацию из mcmc объектов (см. coda), но поиск «принятия» в pdf не дает никаких результатов.


r mcmc
person tomka    schedule 28.09.2017    source источник

Ответы (2)


arrow_upward
1
arrow_downward

Если вам не нравится поведение capture.output, вы можете вместо него использовать sink, который по-прежнему возвращает результаты, но перенаправляет вывод консоли в файл.

sink(file="test.txt")
m <- MCMClogit(y ~ x, burnin = 5000, mcmc = 1000,
               tune = 1.3, B0 = 0, verbose = TRUE)
sink()

out <- readLines("test.txt")
grep( "acceptance rate for beta was", out, value=T)
# [1] "The Metropolis acceptance rate for beta was 0.45533"

class(m)
# [1] "mcmc"

Нет необходимости в assign или <<-.

Если вы хотите получить числовое значение, вы можете извлечь его следующим образом.

ar <- grep( "acceptance rate for beta was", out, value=T)
ar <- as.numeric( gsub("^.* was (.*)", "\\1", ar) )
ar
# [1] 0.45533
person SimonG    schedule 28.09.2017
comment
Если я правильно понимаю, это не сохранит значение коэффициента приема в вектор, верно? - person tomka; 28.09.2017

arrow_upward
1
arrow_downward

Это хак... но, тем не менее, решение:

output = capture.output(MCMClogit(y ~ x, burnin = 5000, mcmc = 1000,
                                  tune = 1.3, B0 = 0, verbose = TRUE))

library(stringr)
library(dplyr)

output %>%
  paste0(collapse = "") %>%
  str_extract("\\d+[.]\\d+(?=[@])") %>%
  as.numeric()

# [1] "0.45533"

Это использует capture.output для сохранения того, что напечатано в консоли из вашего MCMClogit в переменной в виде строки, а затем использует регулярное выражение для извлечения скорости принятия. Что делает регулярное выражение относительно простым, так это то, что скорость принятия окружена цифрами @.

ОП хорошо заметил, что при использовании этого метода MCMClogit придется запускать дважды, что нежелательно, если запуск модели занимает очень много времени. Единственное, что можно сделать, это использовать <<- (может быть опасно) и назначить модель m в глобальной среде:

output = capture.output(m <<- MCMClogit(y ~ x, burnin = 5000, mcmc = 1000,
                                        tune = 1.3, B0 = 0, verbose = TRUE))

Таким образом, объект модели будет сохранен в m, а вывод консоли будет сохранен в output одновременно.

Также обратите внимание, что <<- фактически назначает переменные в родительской среде. В этом случае, когда имеется только одна функция, родительской средой является глобальная среда. Однако в случае, когда есть вложенные функции, родительская среда будет на один уровень выше вложенной, поэтому вместо этого следует использовать assign:

output = capture.output(assign("m", MCMClogit(y ~ x, burnin = 5000, mcmc = 1000,
                                              tune = 1.3, B0 = 0, verbose = TRUE),
                               envir = globalenv()))
person acylam    schedule 28.09.2017
comment
Спасибо - значит, объект mcmc здесь потерян, если я прав? Тогда код нужно будет запускать дважды. - person tomka; 28.09.2017
comment
@tomka Хороший вопрос. Смотрите мои правки. Я бы очень хотел, чтобы кто-нибудь придумал лучший способ, поскольку я часто не решаюсь использовать <<- внутри функции. - person acylam; 28.09.2017
comment
Уточните, чем это опасно? - person tomka; 28.09.2017
comment
@tomka В этом случае все должно быть в порядке, но в некоторых других случаях вы можете неожиданно изменить глобальную переменную, когда вы этого не хотите, или оказаться в ситуации, когда вы ссылаетесь на глобальную переменную внутри функции, когда вы действительно хотели ссылаться на переменную, переданную в качестве параметра. - person acylam; 28.09.2017