如何将 mcmc.list 转换为 bugs 对象？

Posted 2023-02-19

【中文标题】如何将 mcmc.list 转换为 bugs 对象？

【英文标题】：How can I convert an mcmc.list to a bugs object?

【发布时间】：2012-08-18 04:09:43

【问题描述】：

我正在使用rjags R 库。函数coda.samples 产生一个mcmc.list，例如（来自example(coda.samples)）：

library(rjags)
data(LINE)
LINE$recompile()
LINE.out <- coda.samples(LINE, c("alpha","beta","sigma"), n.iter=1000)
class(LINE.out)
[1] "mcmc.list"

但是，我想使用 plot.bugs 函数，它需要 bugs 对象作为输入。

是否可以将对象从mcmc.list 转换为bugs 对象，以便plot.bugs(LINE.out)？

请注意，有一个 similar question on stats.SE 一个多月没有得到答复。该问题已于 2012 年 8 月 29 日结束。

更多提示：

我发现 R2WinBUGS 包有一个函数“as.bugs.array”函数 - 但不清楚该函数如何应用于 mcmc.list。

【问题讨论】：

Abe 为您在 Cross Validated 上的问题提供的答案有什么问题？您能否发布一个图表来显示您想要的上述示例的情节？您在 Cross Validated 上发布了一个图，但它似乎不适用于上面的示例。

@MarkMiller Cross Validated 的答案不完整。

你想要什么具体的加法结果？安倍的答案在我的电脑上运行。了解您想要的附加输出将帮助人们提供必要的代码。这就是为什么我建议您为上面的示例提供一个图表，以准确显示您想要的内容。

在您的 Cross Validated 帖子中，您提供了一个图表，显示了 80% interval for each chain、R-hat 和 medians and 80% intervals 的图表。这就是安倍的回答与您上面的示例所提供的。我在上面的代码中添加的所有内容都是library(R2WinBUGS)，并且我在 Abe 的 plot 声明中添加了一个缺失的括号（我现在已通过提交的编辑将其添加到他的帖子中）。

您在 Cross Validated 上的图包括附加参数图，这可能是因为该图来自不同的示例，或者可能是因为它来自使用相同数据集的不同模型，并且比示例中监控的参数更多更多。这就是为什么我问你想要什么额外的结果。

【参考方案1】：

我不知道这是否会给你你想要的。请注意，model 代码来自使用您的代码，然后在光标处键入 LINE。其余的只是标准错误代码，除了我使用tau = rgamma(1,1) 作为初始值并且不知道它有多标准。我不止一次看到tau = 1 被用作初始值。或许这样会更好。

实际上，我使用您使用的相同model 代码创建了一个rjags 对象，并添加了一个jags 语句来运行它。我承认这与将 coda 输出转换为 bugs 对象不同，但它可能会导致您获得所需的 plot。

如果您只有mcmc.list 而没有model 代码，而您只是想绘制mcmc.list，那么我的回答将无济于事。

library(R2jags)

x <- c(1, 2, 2, 4, 4,  5,  5,  6,  6,  8) 
Y <- c(7, 8, 7, 8, 9, 11, 10, 13, 14, 13) 

N <- length(x)
xbar <- mean(x)

summary(lm(Y ~ x))

x2 <- x - xbar

summary(lm(Y ~ x2))

# Specify model in BUGS language

sink("model1.txt")

cat("

model  
                for( i in 1 : N ) 
                        Y[i] ~ dnorm(mu[i],tau)
                        mu[i] <- alpha + beta * (x[i] - xbar)
                
                tau ~ dgamma(0.001,0.001) 
                sigma <- 1 / sqrt(tau)
                alpha ~ dnorm(0.0,1.0E-6)
                beta ~ dnorm(0.0,1.0E-6)        
        

",fill=TRUE)
sink()

win.data <- list(Y=Y, x=x, N=N, xbar=xbar)

# Initial values
inits <- function() list(alpha=rnorm(1), beta=rnorm(1), tau = rgamma(1,1))

# Parameters monitored
params <- c("alpha", "beta", "sigma")

# MCMC settings
ni <- 25000
nt <-     5
nb <-  5000
nc <-     3

out1 <- jags(win.data, inits, params, "model1.txt", n.chains = nc, 
             n.thin = nt, n.iter = ni, n.burnin = nb)

print(out1, dig = 2)
plot(out1)

#library(R2WinBUGS)
#plot(out1)

编辑：

基于 cmets，也许这样的东西会有所帮助。 str(new.data) 行表明有大量数据可用。如果您只是尝试创建默认图的变体，那么这样做可能只是根据需要提取和子集数据的问题。这里plot(new.data$sims.list$P1) 只是一个简单的例子。在不确切知道您想要什么情节的情况下，我不会尝试更具体的数据提取。如果您发布一个显示您想要的确切情节示例的图形，也许有人可以从这里获取它并发布创建它所需的代码。

顺便说一句，我建议将示例数据集的大小减少到可能是三个链，并且可能不超过 30 次迭代，直到您为您想要的确切绘图获得所需的确切代码：

load("C:/Users/mmiller21/simple R programs/test.mcmc.list.Rdata")

class(test.mcmc.list)

library(R2WinBUGS)

plot(as.bugs.array(sims.array = as.array(test.mcmc.list)))

new.data <- as.bugs.array(sims.array = as.array(test.mcmc.list))

str(new.data)

plot(new.data$sims.list$P1)

编辑：

还要注意：

class(new.data)
[1] "bugs"

而：

class(test.mcmc.list)
[1] "mcmc.list"

这是您帖子的标题所要求的。

【讨论】：

谢谢。与其说我只有一个 mcmc.list，不如说它嵌入在需要 mcmc.list 进行下游处理的代码中（从line 101 here开始。

【参考方案2】：

这不是您问题的解决方案，但作为对您对@andybega 答案的评论的回应，这里有一种将mcmc.list 对象转换为典型的尾文本文件的方法。

mcmc.list.to.coda <- function(x, outdir=tempdir()) 
  # x is an mcmc.list object
  x <- as.array(x)
  lapply(seq_len(dim(x)[3]), function(i) 
    write.table(cbind(rep(seq_len(nrow(x[,,i])), ncol(x)), c(x[,,i])), 
                paste0(outdir, '/coda', i, '.txt'),
                row.names=FALSE, col.names=FALSE)
  )

  cat(paste(colnames(x), 
            1 + (seq_len(ncol(x)) - 1) * nrow(x),
            nrow(x) * seq_len(ncol(x))), 
      sep='\n', 
      file=file.path(outdir, 'codaIndex.txt'))


# For example, using the LINE.out from your question:
mcmc.list.to.coda(LINE.out, tempdir())

【参考方案3】：

不是答案，但这个 blog post 具有以下包装函数，用于使用 R2WinBUGS:::bugs.sims 将尾声输出 (.txt) 转换为 BUGS：

coda2bugs <- function(path, para, n.chains=3, n.iter=5000, 
                      n.burnin=1000, n.thin=2)    
 setwd(path)   
 library(R2WinBUGS)   
 fit <- R2WinBUGS:::bugs.sims(para, n.chains=n.chains, 
        n.iter=n.iter, n.burnin=n.burnin, n.thin=n.thin, 
        DIC = FALSE)   
 class(fit) <- "bugs"   
 return(fit) 
 

【讨论】：

如果我能弄清楚如何将 mcmc.list 保存为 coda 文件，这可能会很有用

参见my answer below 了解将mcmc.list 对象转换为尾声的一种方法。但是，我相信@andybega（Yu-Sung）的解决方案涉及重新拟合和采样。