R中gelman.plot MCMC诊断的in图等价
我正在创建一系列的MCMC诊断图，在r使用am图。我意识到gg中已经有一个软件包可用于MCMC绘图，但其中大部分是为了我自己的教育和实际用途。有一件事我似乎搞不懂，那就是如何在ggplot框架中生成gelman.plot。
gelman.diag函数只返回一个简单的数据点，我想重新创建完整的运行图表，如gelman.plot所示。
是否有人熟悉gelman潜在的降尺度因子的算法结构和/或将其输出移植到ggplot的一种方法？
谢谢!
回答 1
Stack Overflow用户
发布于 2017-01-14 22:43:48
您还没有提供一个可重复的示例，所以我使用了这里的例子。我们需要那个例子中的名为combinedchains的对象。为了避免混乱的答案，我把这方面的代码放在这篇文章的末尾。
现在我们可以在gelman.plot上运行combined.chains了。这就是我们想要复制的情节：
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示例代码
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的形式返回适当的诊断数据。
https://stackoverflow.com/questions/41657817
复制
腾讯云开发者
