gpt4 book ai didi

r - 如何使用plot.mcmc向后验密度图添加垂直线?

转载 作者:行者123 更新时间:2023-12-02 03:15:10 24 4
gpt4 key购买 nike

我经常在已知参数的模拟数据上运行 JAGS 模型。我喜欢 mcmc 对象的默认绘图方法。但是,我想为每个建模参数添加一个 abline(v=TRUE_VALUE) 。这可以让我快速检查后验是否合理。

当然,我可以手动完成此操作,或者可能重新发明轮子并编写我自己的函数。但我想知道是否有一种基于现有 plot 方法的优雅方法。

这是一个有效的示例:

require(rjags)
require(coda)

# simulatee data
set.seed(4444)
N <- 100
Mu <- 100
Sigma <- 15
y <- rnorm(n=N, mean=Mu, sd=Sigma)
jagsdata <- list(y=y)

jags.script <- "
model {
for (i in 1:length(y)) {
y[i] ~ dnorm(mu, tau)
}
mu ~ dnorm(0, 0.001)
sigma ~ dunif(0, 1000)
tau <- 1/sigma^2
}"


mod1 <- jags.model(textConnection(jags.script), data=jagsdata, n.chains=4,
n.adapt=1000)
update(mod1, 200) # burn in
mod1.samples <- coda.samples(model=mod1,
variable.names=c('mu', 'sigma'),
n.iter=1000)
plot(mod1.samples)

enter image description here

我只想为 mu 运行 abline(v=100) ,为 sigma 运行 abline(v=15) 。当然,在许多其他示例中,我会有 5 个、10 个、20 个或更多感兴趣的参数。因此,我对能够为命名参数提供真实值向量感兴趣。

我已经查看了getAnywhere(plot.mcmc)。修改它是一个好方法吗?

最佳答案

好的。所以我修改了 plot.mcmc 如下所示:

my.plot.mcmc <- function (x, trace = TRUE, density = TRUE, smooth = FALSE, bwf, 
auto.layout = TRUE, ask = FALSE, parameters, ...)
{
oldpar <- NULL
on.exit(par(oldpar))
if (auto.layout) {
mfrow <- coda:::set.mfrow(Nchains = nchain(x), Nparms = nvar(x),
nplots = trace + density)
oldpar <- par(mfrow = mfrow)
}
for (i in 1:nvar(x)) {
y <- mcmc(as.matrix(x)[, i, drop = FALSE], start(x),
end(x), thin(x))
if (trace)
traceplot(y, smooth = smooth, ...)
if (density) {
if (missing(bwf)) {
densplot(y, ...); abline(v=parameters[i])
} else densplot(y, bwf = bwf, ...)
}
if (i == 1)
oldpar <- c(oldpar, par(ask = ask))
}
}

然后运行命令

my.plot.mcmc(mod1.samples, parameters=c(Mu, Sigma))

产生这个

enter image description here

请注意,parameters 必须是值的向量,其排序顺序与 JAGS 对变量的排序顺序相同,对于向量来说,这似乎是按字母顺序排列,然后按数字顺序排列。

经验教训

  • 默认情况下,简单地编写一个新的 plot.mcmc 是行不通的,可能是由于命名空间的原因。所以我刚刚创建了一个新函数
  • 我不得不将 set.mfrow 更改为 coda:::set.mfrow 大概也是因为命名空间的原因。
  • 我将 Ask 更改为 ask=FALSE,因为 RStudio 允许浏览数字。

我很高兴听到任何关于覆盖或调整现有 S3 方法的更好方法的建议。

关于r - 如何使用plot.mcmc向后验密度图添加垂直线?,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/10925944/

24 4 0
Copyright 2021 - 2024 cfsdn All Rights Reserved 蜀ICP备2022000587号
广告合作:1813099741@qq.com 6ren.com