1

我需要做两件事。首先，我希望能够在codamcmc 对象中创建从现有变量计算的新变量，以便我可以对新变量运行链诊断。其次，我希望能够在某些 coda plot 函数中索引单个变量，同时仍然查看所有链。

玩具数据。JAGS使用和对睡眠数据进行贝叶斯 t 检验rjags。

data(sleep)

# read in data
y <- sleep$extra
x <- as.numeric(as.factor(sleep$group))
nTotal <- length(y)
nGroup <- length(unique(x)) 
mY <- mean(y)
sdY <- sd(y)

# make dataList
dataList <- list(y = y, x = x, nTotal = nTotal, nGroup = nGroup, mY = mY, sdY = sdY)

# model string
modelString <- "
model{

    for (oIdx in 1:nTotal) {
    y[oIdx] ~ dnorm(mu[x[oIdx]], 1/sigma[x[oIdx]]^2)
    }

    for (gIdx in 1:nGroup) {
    mu[gIdx] ~ dnorm(mY, 1/sdY)
    sigma[gIdx] ~ dunif(sdY/10, sdY*10)
    }
}
" 
writeLines(modelString, con = "tempModel.txt")

# chains

  # 1. adapt
  jagsModel <- jags.model(file = "tempModel.txt",
                          data = dataList,
                          n.chains = 3,
                          n.adapt = 1000)

  # 2. burn-in
  update(jagsModel, n.iter = 1000)

  # 3. generate
  codaSamples <- coda.samples(model = jagsModel,
                              variable.names = c("mu", "sigma"),
                              thin = 15,
                              n.iter = 10000*15/3)

问题一

如果我将coda对象转换为数据框，我可以计算两组估计值之间的差异并绘制这个新变量，就像这样......

df <- as.data.frame(as.matrix(codaSamples))
names(df) <- gsub("\\[|\\]", "", names(df), perl = T) # remove brackets
df$diff <- df$mu1 - df$mu2
ggplot(df, aes(x = diff)) + 
       geom_histogram(bins = 100, fill = "skyblue") + 
       geom_vline(xintercept = mean(df$diff), colour = "red", size = 1, linetype = "dashed")

...但是我如何获得跟踪图？我可以像这样为 coda 对象中的现有变量获取一个...

traceplot(codaSamples[[1]][,1])

...但我希望能够为新diff变量获取它们。

问题二

这让我想到了第二个问题。我希望能够获得单个变量的跟踪图（除其他外）。如上所示，如果我只想查看一个链但我想查看所有链，我可以为单个变量获取它们。我可以使用简单的方法查看模型中所有变量的所有链

情节（codaSamples）

...但是如果我不想或不需要查看所有变量怎么办？如果我只想查看一个甚至两个变量（但不是所有变量）的轨迹图和/或密度图，但图中的所有链都包含在内，该怎么办？

4

0 回答 0