【问题标题】:JAGS error for MCMC Bayesian inferenceMCMC 贝叶斯推理的 JAGS 错误
【发布时间】:2015-12-06 03:25:52
【问题描述】:

在 R 中,我正在对来自 Gamma 分布混合的数据进行 MCMC 贝叶斯推理。这里使用 JAGS。模型文件gmd.bug如下

model {
for (i in 1:N) {
y[i] ~ dsum(p*one, (1-p)*two)
}
one ~ dgamma(alpha1, beta1)
two ~ dgamma(alpha2, beta2)    alpha1 ~ dunif(0, 10)
beta1 ~ dunif(0, 10)
alpha2 ~ dunif(0, 10)
beta2 ~ dunif(0, 10)
p ~ dunif(0, 1)
}

那么,这就是推理阶段

gmd.jags = jags.model("gmd.bug",
data = list(N = NROW(a), y=unlist(a)),
inits = inits, n.chains = 3, n.adapt = 1000)

这是让我困惑的错误

Error in jags.model("gmd.bug", data = list(N = NROW(a), y = unlist(a)),  : 
Error in node y[1]
Node inconsistent with parents

有人知道这里需要修复什么吗?

【问题讨论】:

    标签: r mcmc jags


    【解决方案1】:

    回答 OP 的原始问题y[i] ~ dsum(p*dgamma(alpha1, beta1), (1-p)*dgamma(alpha2, beta2))时,dgamma(alpha1, beta1)需要被[i]索引,如

    gamma1[i] ~ dgamma(alpha1, beta1)
    gamma2[i] ~ dgamma(alpha2, beta2)
    

    回答 OP 的第二个问题(修改后)

    这是你问题的症结所在。但修复它会带来额外的困难,因为为了确保 y[i] 在初始化时与其父级一致,您需要确保在初始化时严格正确 y[i] == p*gamma1[i]+(1-p)*gamma2[i]。如果您让 JAGS 自动处理初始化,它将从先验初始化,而不了解 dsum 对初始值施加的约束,您将收到错误消息。一种策略是在y 初始化gamma1gamma2

    以下代码适用于我(但当然你会想要运行更多的迭代):

    # Data simulation:
    library(rjags)
    N=200
    alpha1 <- 3
    beta1 <- 3
    alpha2 <- 5
    beta2 <- 1
    p <- .7
    
    y <- vector(mode="numeric", length=N)
    for(i in 1:N){
      y[i] <- p*rgamma(1,alpha1,beta1) + (1-p)*rgamma(1,alpha1,beta1)
    }
    
    # JAGS model
    sink("mymodel.txt")
    cat("model{
        for (i in 1:N) {
          gamma1[i] ~ dgamma(alpha1, beta1)
          gamma2[i] ~ dgamma(alpha2, beta2)
          pg1[i] <- p*gamma1[i]
          pg2[i] <- (1-p)*gamma2[i]
          y[i] ~ dsum(pg1[i], pg2[i])
        }
        alpha1 ~ dunif(0, 10)
        beta1 ~ dunif(0, 10)
        alpha2 ~ dunif(0, 10)
        beta2 ~ dunif(0, 10)
        p ~ dunif(0, 1)
        }", fill=TRUE)
    sink()
    jags.data <- list(N=N, y=y)
    
    inits <- function(){list(gamma1=y, gamma2=y)}
    
    params <- c("alpha1", "beta1", "alpha2", "beta2", "p")
    
    nc <- 5
    n.adapt <-200
    n.burn <- 200
    n.iter <- 1000
    thin <- 10
    mymodel <- jags.model('mymodel.txt', data = jags.data, inits=inits, n.chains=nc, n.adapt=n.adapt)
    update(mymodel, n.burn)
    mymodel_samples <- coda.samples(mymodel,params,n.iter=n.iter, thin=thin)
    

    【讨论】:

    • 我刚刚发现存在一个问题,即 gamma1 和 gamma2 在后验中几乎相同,并且 p 的估计值没有收敛。我尝试了几种情况,这似乎是一个问题。
    • 如果 p 不能收敛,那么 gamma1 和 gamma2 是相同的是有道理的(因为如果模型无法区分 p=Pp=1-P 那么显然将无法学习 gamma1 和gamma2 不同)。有两个可能的问题来源。一是混合很差,在这种情况下,您可以尝试运行模型更长时间,看看 p 是否最终收敛。另一个是您没有足够的数据来估计模型,在这种情况下,您可以尝试模拟一个更大的数据集(y 中的元素更多),看看模型是否可以用更多的数据进行估计。
    • 谢谢。但是这里的 y 对我来说是一个给定的数据集(200 obs,还不够),我运行了 10000 次迭代。看到五个 MCMC 链都在暗示不同的 p,这很有趣。不管怎样,我从你身上学到了很多。
    • 啊!如果所有链都建议不同的 p,这表明问题出在混合上,而不是模型的可估计性!尝试运行(更多)更长的时间,或者您可以考虑在 Stan 而不是 JAGS 中拟合模型。祝你好运!
    • 考虑到您的模型的一些额外细节,您可以选择能够真正为其提供良好开端的初始值(这没有问题,因为初始值应该独立于 MCMC 最终提出的任何位置)
    猜你喜欢
    • 1970-01-01
    • 2015-10-09
    • 1970-01-01
    • 1970-01-01
    • 2013-05-28
    • 1970-01-01
    • 2018-12-01
    • 2017-11-14
    • 1970-01-01
    相关资源
    最近更新 更多