如何将mcmc.list转换为bug对象? [英] How can I convert an mcmc.list to a bugs object?

问题描述

我正在使用rjags R库.函数coda.samples生成一个mcmc.list，例如(来自example(coda.samples)):

I am using the rjags R library. The function coda.samples produces an mcmc.list, for example (from 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对象作为输入.

However, I would like to use the plot.bugs function, which requires a bugs object as input.

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

Is it possible to convert an object from an mcmc.list to a bugs object, so that plot.bugs(LINE.out)?

请注意，有关stats.SE的类似问题，已有一个多月没有得到解答.这个问题的悬赏金在2012年8月29日结束.

Note that there is a similar question on stats.SE that has been unanswered for over a month. That question had a bounty that ended on 08/29/2012.

我发现R2WinBUGS软件包具有函数"as.bugs.array"功能-但目前尚不清楚该功能如何应用于mcmc.list.

I have discovered that the R2WinBUGS package has a function "as.bugs.array" function - but it is not clear how the function can be applied to an mcmc.list.

推荐答案

我不知道这是否可以满足您的需求.请注意，model代码来自于使用您的代码，然后在光标处键入LINE.其余的只是标准的bug代码，除了我使用tau = rgamma(1,1)作为初始值，并且不知道它的标准是什么.我不止一次看到tau = 1用作初始值.也许那会更好.

I do not know whether this will give you what you want. Note that the model code came from using your code and then typing LINE at the cursor. The rest is just standard bugs code, except I used tau = rgamma(1,1) for an initial value and do not know how standard that is. More than once I have seen tau = 1 used as an initial value. Perhaps that would be better.

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

In effect, I created an rjags object using the same model code you were using and added a jags statement to run it. I admit that is not the same thing as converting coda output to a bugs object, but it might result in you getting the desired plot.

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

If all you have is an mcmc.list and no model code and you simply want to plot the mcmc.list, then my answer will not help.

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)只是一个简单的示例.在不确切知道您想要什么绘图的情况下，我将不会尝试更具体的数据提取.如果您发布的图显示了确切类型的绘图示例，则您可能希望有人可以从此处获取它并发布创建它所需的代码.

Based on the comments perhaps something like this will help. The line str(new.data) suggests that a large amount of data are available. If you are simply trying to create variations of default plots then doing so may only be a matter of extracting and subsetting the data as desired. Here plot(new.data$sims.list$P1) is just one straight-forward example. Without knowing exactly what plot you want I will not attempt more specific data extractions. If you post a figure showing an example of the exact kind of plot you want perhaps someone can take it from here and post the code needed to create it.

顺便说一句，我建议将示例数据集的大小减少到大约三个链，并且不超过30次迭代，直到您拥有所需的精确绘图所需的确切代码为止:

By the way, I recommend reducing the size of the example data set to perhaps three chains and perhaps no more than 30 iterations until you have the exact code you want for the exact plot you want:

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"

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

which is what the title of your post requests.

这篇关于如何将mcmc.list转换为bug对象?的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！