ggplot 相当于 gelman.plot MCMC 诊断在 r
ggplot equivalent of gelman.plot MCMC diagnostic in r
我正在使用 ggplot 在 r 中创建一系列 MCMC 诊断图。我意识到 gg 中已经有一个用于 MCMC 绘图的包,但其中大部分是为了我自己的教育和实际使用。我似乎无法弄清楚的一件事是如何在 ggplot 框架中生成 gelman.plot。
gelman.diag 函数只有 returns 一个简单的数据点,我想重新创建完整的 运行 图表,如 gelman.plot 所示。
是否有人熟悉 gelman 潜在比例缩减因子的算法结构 and/or 一种将其输出移植到 ggplot 的方法?
谢谢!
您没有提供可重现的示例,所以我使用了 example here。我们需要该示例中名为 combinedchains
的对象。为了避免混淆答案,我把代码放在了这个 post.
的末尾
现在我们可以在 combined.chains
上 运行 gelman.plot
。这是我们要复制的情节:
library(coda)
gelman.plot(combined.chains)
要创建 ggplot 版本,我们需要获取绘图数据。我以前没有做过 MCMC,所以我打算让 gelman.plot
为我生成数据。对于您的实际用例,您可能可以直接生成适当的数据。
让我们看看 gelman.plot
在做什么:我们可以通过在控制台中键入裸函数名称来查看该函数的代码。部分功能代码如下。 ...
显示我为简洁起见删除了原始代码部分的位置。请注意对 gelman.preplot
的调用,该函数的输出存储在 y
中。另请注意,y
在末尾被无形地 return 编辑。 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
return不可见的数据并将其存储在一个对象中:
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)
。然后你可以看到这个函数在做什么,如果你愿意,可以构建你自己的函数,以 return 更适合 ggplot 的形式的适当诊断数据。
我正在使用 ggplot 在 r 中创建一系列 MCMC 诊断图。我意识到 gg 中已经有一个用于 MCMC 绘图的包,但其中大部分是为了我自己的教育和实际使用。我似乎无法弄清楚的一件事是如何在 ggplot 框架中生成 gelman.plot。
gelman.diag 函数只有 returns 一个简单的数据点,我想重新创建完整的 运行 图表,如 gelman.plot 所示。
是否有人熟悉 gelman 潜在比例缩减因子的算法结构 and/or 一种将其输出移植到 ggplot 的方法?
谢谢!
您没有提供可重现的示例,所以我使用了 example here。我们需要该示例中名为 combinedchains
的对象。为了避免混淆答案,我把代码放在了这个 post.
现在我们可以在 combined.chains
上 运行 gelman.plot
。这是我们要复制的情节:
library(coda)
gelman.plot(combined.chains)
要创建 ggplot 版本,我们需要获取绘图数据。我以前没有做过 MCMC,所以我打算让 gelman.plot
为我生成数据。对于您的实际用例,您可能可以直接生成适当的数据。
让我们看看 gelman.plot
在做什么:我们可以通过在控制台中键入裸函数名称来查看该函数的代码。部分功能代码如下。 ...
显示我为简洁起见删除了原始代码部分的位置。请注意对 gelman.preplot
的调用,该函数的输出存储在 y
中。另请注意,y
在末尾被无形地 return 编辑。 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
return不可见的数据并将其存储在一个对象中:
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)
。然后你可以看到这个函数在做什么,如果你愿意,可以构建你自己的函数,以 return 更适合 ggplot 的形式的适当诊断数据。