Michael Platzer
2014-03-04
Hello,
I am a new JAGS user, and got unfortunately quickly stuck in a problem which I can't seem to solve, and hope for your help now.
Here is a (reduced) description of my problem: Data is generated in two steps. first a draw from a poisson distribution x ~ Pois(lambda)
, and then a draw from a gamma distribution, whereas x
is used as the shape parameter.
When trying to re-identify the parameters with JAGS, an error is thrown which seems to indicate that JAGS cannot handle 0 values for the shape parameter (see code A below). I then tried to implement a workaround with an if/else construct as suggested here http://stackoverflow.com/questions/15414303/choosing-different-distributions-based-on-if-else-condition-in-winbugs-jags but failed as well (see code B below).
Can you please confirm, whether it is indeed true that the gamma distribution cannot handle shape 0, and if so, would you happen to know of a workaround in this case?
many thanks,
Michael
code A:
# generate artificial sample set.seed(1) N <- 1000 x <- rpois(N, lambda=2) t <- rgamma(N, x, rate=4) model <- function() { for (i in 1:N) { x[i] ~ dpois(lambda) t[i] ~ dgamma(x[i], rate) } lambda ~ dunif(1, 5) rate ~ dunif(1, 5) } R2jags::jags(data=c("x", "t", "N"), inits=function() list(lambda=3, rate=3), c("lambda", "rate"), model.file=model) # module glm loaded # Compiling model graph # Resolving undeclared variables # Allocating nodes # Graph Size: 2005 # # Initializing model # Deleting model # # Error in jags.model(model.file, data = data, inits = init.values, n.chains = n.chains, : # Error in node t[10] # Invalid parent values
code B: workaround attempt
model <- function() { for (i in 1:N) { x[i] ~ dpois(lambda) tt[i,1] <- 0 tt[i,2] ~ dgamma(x[i], rate) tif[i] <- 1 + step(x[i]-1) t[i] ~ tt[i, tif[i]] } lambda ~ dunif(1, 5) rate ~ dunif(1, 5) } R2jags::jags(data=c("x", "t", "N"), inits=function() list(lambda=3, rate=3), c("lambda", "rate"), model.file=model) # Error in jags.model(model.file, data = data, inits = init.values, n.chains = n.chains, : # Error parsing model file: # syntax error on line 8 near "tt"
Martyn Plummer
2014-03-04
There is no such thing as a gamma distribution with zero shape parameter. It simply doesn't exist.
Michael Platzer
2014-03-04
Valid point. And I guess it is wrong to expect from JAGS to mimic the non-standard behavior of R in that respect.
Still, the problem then remains of how to define different distributions depending on the value of x
. See my attempt in code snippet B.
The background is, that I want to model the elapsed time t
until x
events have occurred, whereas waiting times between events are exponentially distributed. t
thus follows an Erlang-x
distribution, but is for the edge case of x=0
zero.
Martyn Plummer
2014-03-04
These two-stage models do require degenerate distributions to be defined in JAGS. At the request of users, I have added degenerate cases for the binomial, negative binomial and Poisson distributions. But I'm really reluctant to do it for a continuous distribution because it creates an infinite density, e.g. in R you have
> dgamma(0, shape=0, rate=1) [1] Inf > dgamma(0, shape=0, rate=10) [1] Inf > dgamma(0, shape=0, rate=100) [1] Inf
If you are trying sample the rate parameter, then the presence of an infinite likelihood is rather awkward.
You can work around the problem by adding a small quantity to the shape parameter to keep it away from zero:
t[i] ~ dgamma(x[i] + 1.0E-3, rate)
Then you must set t[i]
to missing whenever x[i]==0
to ensure that t[i]
is uninformative about rate
in such cases.
t <- rgamma(N, x, rate=4) t[x==0] <- NA
Ultimately, it would be useful to have if/else constructs in the BUGS language that work whenever the predicate is a function of the data and can be evaluated at compile time.
Then you could have
if (x[i] > 0) {
t[i] ~ dgamma(x[i], rate)
}
This is on the list of things to do when I overhaul the compiler.
Michael Platzer
2014-03-05
Excellent, this helps a lot!
I actually already tried out setting t[x==0]<-NA
, and and also incrementing x
by a small quantity, but I did so separately. Only when doing both, the jags simulation compiles and runs properly.