Übersetzen WinBUGS Modell JAGS (mit R)

Question

Ich versuche, das folgende Modell für WinBUGS in JAGS geschrieben umzusetzen:

model { 
    for (i in 1:N) {       
    wtp[i] ~ dweib(r[G[i]], mu[i])I(lower[i], upper[i]) 
    mu[i] <- exp(beta[G[i]]) 
    G[i] ~ dcat(P[])  
    }          
    P[1] ~ dunif(0.01, 0.99) 
    P[2] <- 1 - P[1] 
    r[1] ~ dunif(1, 10) 
    r[2] ~ dunif(0.1, 10) 
    beta[1] ~ dunif(0, 1000)  
    beta[2] ~ dunif(-1000, 0)     
    weibmed[1] <- pow(log(2) * exp(-beta[1]), 1/r[1]) 
    weibmed[2] <- pow(log(2) * exp(-beta[2]), 1/r[2]) 
    weibmed[3] <- pow(log(1/(1 - 0.5 + P[1])) * exp(-beta[2]), 1/r[2]) 
    weibmean[1] <- pow(exp(-beta[1]), 1/r[1]) * exp(loggam((1 + r[1])/r[1])) 
    weibmean[2] <- pow(exp(-beta[2]), 1/r[2]) * exp(loggam((1 + r[2])/r[2])) 
    weibmean[3] <- P[1] * weibmean[1] + P[2] * weibmean[2] 
} 

Ich dachte, es einfach wäre, es in JAGS in Gang zu bringen mit:

library(rjags) 

txt <- 'model { 
    for (i in 1:N) {       
    wtp[i] ~ dweib(r[G[i]], mu[i])T(lower[i], upper[i]) 
    mu[i] <- exp(beta[G[i]]) 
    G[i] ~ dcat(P[])  
    }          
    P[1] ~ dunif(0.01, 0.99) 
    P[2] <- 1 - P[1] 
    r[1] ~ dunif(1, 10) 
    r[2] ~ dunif(0.1, 10) 
    beta[1] ~ dunif(0, 1000)  
    beta[2] ~ dunif(-1000, 0)     
    weibmed[1] <- pow(log(2) * exp(-beta[1]), 1/r[1]) 
    weibmed[2] <- pow(log(2) * exp(-beta[2]), 1/r[2]) 
    weibmed[3] <- pow(log(1/(1 - 0.5 + P[1])) * exp(-beta[2]), 1/r[2]) 
    weibmean[1] <- pow(exp(-beta[1]), 1/r[1]) * exp(loggam((1 + r[1])/r[1])) 
    weibmean[2] <- pow(exp(-beta[2]), 1/r[2]) * exp(loggam((1 + r[2])/r[2])) 
    weibmean[3] <- P[1] * weibmean[1] + P[2] * weibmean[2] 
}' 

set.seed(3.14159) 
dat <- list(N = 1000, lower = rep(0, 1000), upper = runif(1000, 5, 200000)) 
ini <- list(P = c(0.4, NA), r = c(8.2, 1.2), beta = c(3.8, -6.5)) 

mod <- jags.model(
    file = textConnection(txt), 
    data = dat, 
    inits = c(ini, .RNG.name = 'base::Mersenne-Twister', .RNG.seed = 314159), 
    n.chains = 1, 
    n.adapt = 100 
) 

sam.jags <- coda.samples(
    model = mod, 
    variable.names = c('P', 'r', 'beta', 'weibmed', 'weibmean'), 
    n.iter = 400, 
    n.thin = 1 
) 

durch nur die I() durch eine T() ersetzen. Dies erzeugt die coda.samples() Fehler:

Error: Error in node weibmed[3] 
Invalid parent values 

Wenn ich die Überwachung von weibmed und weibmean dann coda.samples() Werke ignorieren aber die Parameterschätzungen:

   Mean   SD Naive SE Time-series SE 
P[1]  0.4840704 0.2769491 0.01384746  0.01384746 
P[2]  0.5159296 0.2769491 0.01384746  0.01384746 
beta[1] 509.3614647 295.0860473 14.75430237 14.75430237 
beta[2] -487.5362940 285.4126899 14.27063449 14.27063449 
r[1]  5.2054730 2.6330434 0.13165217  0.13165217 
r[2]  5.0478143 2.9480476 0.14740238  0.14740238 

nicht vergleichbar sind diejenigen, die ich bekomme, wenn WinBUGS mit:

library(R2WinBUGS) 

sam.bugs <- bugs(
    model.file = 'model.bug', 
    data = dat, 
    inits = list(ini), 
    parameters.to.save = c('P', 'r', 'beta'), #, 'weibmed', 'weibmean'), 
    n.chains = 1, 
    n.burnin = 100, 
    n.iter = 500, 
    n.thin = 1, 
    debug = F, 
    DIC = F, 
    bugs.seed = 314159 
) 

Inference for Bugs model at "3mixout2.bug", fit using WinBUGS, 
1 chains, each with 500 iterations (first 100 discarded) 
n.sims = 400 iterations saved 
     mean sd 2.5% 25% 50% 75% 97.5% 
P[1]  0.4 0.0 0.3 0.4 0.4 0.4 0.4 
P[2]  0.6 0.0 0.6 0.6 0.6 0.6 0.7 
r[1]  7.2 0.8 6.2 6.6 6.9 7.7 9.3 
r[2]  1.5 0.0 1.4 1.4 1.5 1.5 1.6 
beta[1] 5.5 0.5 4.6 5.0 5.4 5.8 6.6 
beta[2] -7.2 0.2 -7.5 -7.3 -7.2 -7.1 -6.9 

Irgendwelche Gedanken oder Vorschläge?

Answer 1

JAGS Handbuch sagt;

pow(x, z) || Power function || Real || If x < 0 then z is integer

Wenn 0.5 < P[1], log(1/(1 - 0.5 + P[1])) * exp(-beta[2]) < 0, so weibmed[3], pow(negative, non_integer), NaN wird. Soweit ich weiß, erlaubt coda.samples() keine überwachten Variablen zu nehmen NaN.

Wenn Sie P[1] ~ dunif(0.01, 0.5) oder außer weibmed[3] aus der Liste der überwachten Variablen verwenden, indem Sie den Namen ändern, z. B. weibmed3 <- pow(log(1/(1 - 0.5 + P[1])) * exp(-beta[2]), 1/r[2]), wird Ihr Code ausgeführt.

Answer 2

Es scheint, dass dies ein Fall ist, wo dinterval sollte statt T() verwendet werden:

txt2 <- ' 
data { 
    x <- rep(1, N) 
} 
model { 
    for (i in 1:N) {  
    x[i] ~ dinterval(wtp[i], B[i, ])     
    wtp[i] ~ dweib(r[G[i]], mu[i]) 
    mu[i] <- exp(beta[G[i]]) 
    G[i] ~ dcat(P[]) 
    }          
    P[1] ~ dunif(0, 1) 
    P[2] <- 1 - P[1] 
    r[1] ~ dunif(0, 10) 
    r[2] ~ dunif(0, 10) 
    beta[1] ~ dunif(1, 1000)  
    beta[2] ~ dunif(-1000, 0)     
    weibmed[1] <- pow(log(2) * exp(-beta[1]), 1/r[1]) 
    weibmed[2] <- pow(log(2) * exp(-beta[2]), 1/r[2]) 
    # weibmed[3] <- pow(log(1/(1 - 0.5 + P[1])) * exp(-beta[2]), 1/r[2]) 
    weibmean[1] <- pow(exp(-beta[1]), 1/r[1]) * exp(loggam((1 + r[1])/r[1])) 
    weibmean[2] <- pow(exp(-beta[2]), 1/r[2]) * exp(loggam((1 + r[2])/r[2])) 
    # weibmean[3] <- P[1] * weibmean[1] + P[2] * weibmean[2] 
}' 

mod <- jags.model(
    file = textConnection(txt2), 
    data = list(N = dat[[1]], B = cbind(dat$lower, dat$upper)), 
    inits = c(ini, .RNG.name = 'base::Mersenne-Twister', .RNG.seed = 314159), 
    n.chains = 1, 
    n.adapt = 100 
)