Ошибка JAGS для MCMC Байесовский вывод

Question

В R я использую байесовский вывод MCMC для данных из смеси Gamma-дистрибутивов. Здесь используется 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 

Каждый знает, что должно быть исправлено здесь?

Answer 1

Ответ на первоначальный вопрос OP в Когда вы пишете y[i] ~ dsum(p*dgamma(alpha1, beta1), (1-p)*dgamma(alpha2, beta2)), dgamma(alpha1, beta1) должен быть проиндексированы [я], как и в

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

Ответ на второй вопрос OP (после редактирования)

В этом суть вашей проблемы. Но исправление этого вызывает дополнительные трудности, потому что для обеспечения того, чтобы y [i] был совместим с родителями при инициализации, вам нужно убедиться, что при инициализации строго верно, что y[i] == p*gamma1[i]+(1-p)*gamma2[i]. Если вы позволяете JAGS обрабатывать инициализацию автоматически, она будет инициализироваться из priors, не понимая ограничения на начальные значения, налагаемые dsum, и вы получите сообщение об ошибке. Одна из стратегий заключается в том, чтобы инициализировать как gamma1, так и gamma2 по адресу y.

Следующий код работает для меня (но, конечно, вы хотите запустить много больше итераций):

# 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)