gamma distribution with shape 0

Help
2014-03-04
2014-03-05
  • 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
    Martyn Plummer
    2014-03-04

    There is no such thing as a gamma distribution with zero shape parameter. It simply doesn't exist.

     
  • 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
    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.

     
  • 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.