3

我经常在具有已知参数的模拟数据上运行 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)

在此处输入图像描述

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

我看过了getAnywhere(plot.mcmc)。修改它是一个好方法吗？

4

1 回答 1

7

好的。所以我修改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))

产生这个

在此处输入图像描述

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

得到教训

  • 默认情况下，简单地编写一个 newplot.mcmc不起作用大概是因为命名空间。所以我刚刚创建了一个新功能
  • 我不得不更改set.mfrow为coda:::set.mfrow大概也是因为命名空间。
  • 我将 ask 更改为ask=FALSE，因为 RStudio 允许浏览数字。

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

于 2012-06-07T06:32:48.050 回答