ggplot 相当于 r 中的 gelman.plot MCMC 诊断

Question

我正在使用 ggplot 在 r 中创建一系列 MCMC 诊断图。我意识到 gg 中已经有一个用于 MCMC 绘图的包，但其中大部分是为了我自己的教育和实际使用。我似乎无法弄清楚的一件事是如何在 ggplot 框架中生成 gelman.plot。

gelman.diag 函数只返回一个简单的数据点，我想重新创建完整的运行图表，如 gelman.plot 所示。

是否有人熟悉格尔曼潜在比例缩减因子的算法结构和/或将其输出移植到 ggplot 的方法？

谢谢！

Answer 1

您没有提供可重现的示例，所以我使用了example here。我们需要该示例中名为combinedchains 的对象。为了避免混淆答案，我已将代码放在本文末尾。

现在我们可以在combined.chains 上运行gelman.plot。这是我们要复制的情节：

library(coda)
gelman.plot(combined.chains)

要创建 ggplot 版本，我们需要获取绘图的数据。我之前没有做过 MCMC，所以我打算让gelman.plot 为我生成数据。对于您的实际用例，您可能可以直接生成适当的数据。

让我们看看gelman.plot 在做什么：我们可以通过在控制台中输入裸函数名来查看该函数的代码。部分功能代码如下。为简洁起见，... 显示了我删除原始代码部分的位置。注意对gelman.preplot 的调用，该函数的输出存储在y 中。另请注意，y 在最后不可见地返回。 y 是一个列表，其中包含我们在 ggplot 中创建 gelman.plot 所需的数据。

gelman.plot = function (x, bin.width = 10, max.bins = 50, confidence = 0.95, 
          transform = FALSE, autoburnin = TRUE, auto.layout = TRUE, 
          ask, col = 1:2, lty = 1:2, xlab = "last iteration in chain", 
          ylab = "shrink factor", type = "l", ...) 
{ 
  ...
  y <- gelman.preplot(x, bin.width = bin.width, max.bins = max.bins, 
                      confidence = confidence, transform = transform, autoburnin = autoburnin)
...
  return(invisible(y))
}

那么，让我们获取gelman.plot 隐形返回的数据并将其存储在一个对象中：

gp.dat = gelman.plot(combinedchains)

现在是 ggplot 版本。首先，gp.dat 是一个列表，我们需要将该列表的各个部分转换为 ggplot 可以使用的单个数据框。

library(ggplot2)
library(dplyr)
library(reshape2)

df = data.frame(bind_rows(as.data.frame(gp.dat[["shrink"]][,,1]), 
                          as.data.frame(gp.dat[["shrink"]][,,2])), 
                q=rep(dimnames(gp.dat[["shrink"]])[[3]], each=nrow(gp.dat[["shrink"]][,,1])),
                last.iter=rep(gp.dat[["last.iter"]], length(gp.dat)))

对于情节，我们将df 融合为长格式，以便我们可以将每个链放在一个单独的方面。

ggplot(melt(df, c("q","last.iter"), value.name="shrink_factor"), 
       aes(last.iter, shrink_factor, colour=q, linetype=q)) + 
  geom_hline(yintercept=1, colour="grey30", lwd=0.2) +
  geom_line() +
  facet_wrap(~variable, labeller= labeller(.cols=function(x) gsub("V", "Chain ", x))) +
  labs(x="Last Iteration in Chain", y="Shrink Factor",
       colour="Quantile", linetype="Quantile") +
  scale_linetype_manual(values=c(2,1))

---

用于创建combinedchains 对象的MCMC 示例代码（从here 复制的代码）：

trueA = 5
trueB = 0
trueSd = 10
sampleSize = 31

x = (-(sampleSize-1)/2):((sampleSize-1)/2)
y =  trueA * x + trueB + rnorm(n=sampleSize,mean=0,sd=trueSd)

likelihood = function(param){
  a = param[1]
  b = param[2]
  sd = param[3]

  pred = a*x + b
  singlelikelihoods = dnorm(y, mean = pred, sd = sd, log = T)
  sumll = sum(singlelikelihoods)
  return(sumll)
}

prior = function(param){
  a = param[1]
  b = param[2]
  sd = param[3]
  aprior = dunif(a, min=0, max=10, log = T)
  bprior = dnorm(b, sd = 5, log = T)
  sdprior = dunif(sd, min=0, max=30, log = T)
  return(aprior+bprior+sdprior)
}

proposalfunction = function(param){
  return(rnorm(3,mean = param, sd= c(0.1,0.5,0.3)))
}

run_metropolis_MCMC = function(startvalue, iterations) {
  chain = array(dim = c(iterations+1,3))
  chain[1,] = startvalue
  for (i in 1:iterations) {
    proposal = proposalfunction(chain[i,])

    probab = exp(likelihood(proposal) + prior(proposal) - likelihood(chain[i,]) - prior(chain[i,]))

    if (runif(1)  <  probab){
      chain[i+1,] = proposal
    }else{
      chain[i+1,] = chain[i,]
    }
  }
  return(mcmc(chain))
}

startvalue = c(4,2,8)
chain = run_metropolis_MCMC(startvalue, 10000)

chain2 = run_metropolis_MCMC(startvalue, 10000)
combinedchains = mcmc.list(chain, chain2)

---

更新： gelman.preplot 是一个内部 coda 函数，用户无法直接看到。要获取功能代码，请在控制台中输入getAnywhere(gelman.preplot)。然后您可以查看该函数正在做什么，如果您愿意，可以构建自己的函数以更适合 ggplot 的形式返回适当的诊断数据。