Anonymous
2011-06-22
Attempting to set a monitor on the variable 'phi' in the following model:
data { for (i in 1:N){ for (j in 1:P[ i ]){ zeros[i,j] <- 0 } } } model { #hyperpriors a ~ dnorm(60,0.25)I(0,) b ~ dgamma(1,16)I(0,) q ~ dgamma(1200,10) r ~ dgamma(40,20) #Outer loop for (i in 1:N){ #priors alpha[ i ] ~ dnorm(a,b)I(0,) theta[ i ] ~ dgamma(q,r)I(1.01,) #Inner loop for (j in 1:P[ i ]){ #'Zeros' trick to produce custom likelihood zeros[i,j] ~ dpois(phi[i,j]) #log-likelihood phi[i,j] <- -(log(theta[ i ]) - (theta[ i ]*log(alpha[ i ])) + ((theta[ i ]-1)*log(D[i,j])))*delta[ i ] + (1/P[ i ])*(pow((Age[ i ]+1)/alpha[ i ],theta[ i ])-pow(Age[ i ]/alpha[ i ],theta[ i ])) } } avgAlph <- sum(alpha[])/N avgThet <- sum(theta[])/N } }
Yields the error "Failed to set trace monitor for node phi". Additional
information:
-I cannot post a copy of my data for confidentiality reasons
-There are no undefined values in the data D, nor in the other defined quantities P and Age; Am I wrong to think that phi is therefore defined for all values?
-It doesn't seem to be an indexing problem, as I can happily set monitors for alpha and theta
-The simulation seems to be running just fine in all other respects
My current hypothesis: the 'zeros trick' to declare a custom likelihood
function in WinBUGS has to be handled differently in JAGS. I can't imagine I
am the first person to attempt this, but I haven't found any material
discussing it either.
Martyn Plummer
2011-06-23
Sorry I can't do much without a worked example: model, initial values, data &
script. You can send these by mail and I will treat them as confidential. Fake
data will do: I just need something that will reproduce the problem.
Anonymous
2011-06-30
The email I sent seems to have bounced; here is a dummy script/data that
reproduces the problem using the above model:
Script
library(rjags) foo <- function(file, data, variable.names, n.iter, n.adapt, thin) { m <- jags.model(file, data, n.chains = 1, n.adapt = n.adapt) coda.samples(m, variable.names, n.iter, thin) } ## set data file <- "Test.bug" data <- read.data("TestData.txt","bugs") ## set parameters vars <- c("phi", "theta") n.iter <- 15000 n.adapt <- 3000 thin <- 5 n.chains <- 4 params <- list("foo", "jags.model", "coda.samples", "mcmc", "mcmc.list", "file", "data", "vars", "n.iter", "n.adapt", "thin") post <- mcmc.list(foo(file, data, vars, n.iter, n.adapt, thin))
Data
list(D = structure( .Data=c(99,99, 5.8,99, 99,99, 99,99, 10.6,99, 99,99, 99,99, 3.6,12.1), .Dim=c(8,2) ), delta=c(0,1,0,0,1,0,0,1), Age=c(6,5,12,14,10,3,6,3), P=c(1,1,1,1,1,1,1,2), N=8 )
Martyn Plummer
2011-06-30
You cannot monitor phi because you have not defined every element.
Anonymous
2011-07-04
I think I understand: the problem is I'm iterating to j=2 for some values of
phi, but only j=1 for others, leaving empty entries in the array phi, correct?
Martyn Plummer
2011-07-05
Yes that's right. You can fill in the blanks with any value you wish, e.g.
for (k in 1:N) { for (j in (P[k]+1):PMAX) { phi[k,j] <- 0 } }
where PMAX=2. Then phi will be fully defined and you can monitor it.