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

Question

我正在使用 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。

Answer 1

我不知道这是否会给你你想要的。请注意 model代码来自使用您的代码，然后输入 LINE在光标处。其余的只是标准的错误代码，除了我使用了 tau = rgamma(1,1)对于一个初始值，不知道那是多么标准。我不止一次见过 tau = 1用作初始值。或许那样会更好。

实际上，我创建了一个 rjags对象使用相同的 model您正在使用的代码并添加了 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)

编辑:

根据评论，也许这样的事情会有所帮助。线路 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"

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