有没有办法保存MCMClogit使用的MCMC算法的接受率？

Question

R 包 MCMCpack 通过 MCMClogit 提供贝叶斯逻辑回归。此函数打印 MCMC (Metropolis-Hastings) 算法的接受率，如果 verbose=TRUE 但似乎 return 它在 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) returning NULL 或 names(attributes(m)) 中找不到此信息 returns

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

帮助文件表明 coda 包允许从 mcmc 对象（参见 coda）中汇总信息，但在 pdf 中搜索 'acceptance' 不会产生任何结果结果。

Answer 1

这是一个 hack...但仍然是一个解决方案：

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 在控制台中打印的内容作为字符串存储在变量中，然后使用正则表达式提取接受率。使正则表达式相对简单的原因是接受率被 @ 包围。

OP 提出了一个很好的观点，即通过使用此方法，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()))

Answer 2

如果您不喜欢 capture.output 的行为，您可以改用 sink，它仍然是 returns 结果，但会将控制台输出重定向到一个文件中。

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

有没有办法保存MCMClogit使用的MCMC算法的接受率？

Is there a way to save the acceptance rate of the MCMC algorithm used by MCMClogit?

r

mcmc