Specify a discrete weibull distribution in JAGS or BUGS for R

Question

I am fitting a weibull model to discrete values using JAGS in R. I have no problem fitting a weibull to continuous data, but I run in to trouble when I switch to discrete values.

Here is some data, and code to fit a weibull model in JAGS:

#draw data from a weibull distribution
y <- rweibull(200, shape = 1, scale = 0.9)
#y <- round(y)

#load jags, specify a jags model.
library(runjags)

j.model ="
model{
for (i in 1:N){
y[i] ~ dweib(shape[i], scale[i])

shape[i] <- b1
scale[i] <- b2
}

#priors
b1 ~ dnorm(0, .0001) I(0, )
b2 ~ dnorm(0, .0001) I(0, )
}
"

#load data as list
data <- list(y=y, N = length(y))

#run jags model.
jags.out <-     run.jags(j.model,
                         data=data,
                         n.chains=3,
                         monitor=c('b1','b2')
                        )
summary(jags.out)

This model fits fine. However, if I transform y values to discrete values using y <- round(y), and run the same model, it fails with the error Error in node y[7], Node inconsistent with parents. The particular number of the node changes every time I try, but its always a low number.

I know I can make this run by adding a very small number to all of my values, however, this does not account for the fact that the data are discrete. I know discrete weibull distributions exists, but how can I implement one in JAGS?

Answer 1

You can use the 'ones trick' to implement a discrete weibull distribution in JAGS. Using the pmf here we can make a function to generate some data:

pmf_weib <- function(x, scale, shape){

  exp(-(x/scale)^shape) - exp(-((x+1)/scale)^shape)
}

# probability of getting 0 through 200 with scale = 7 and shape = 4
probs <- pmf_weib(seq(0,200), 7, 4) 

y <- sample(0:200, 100, TRUE, probs ) # sample from those probabilities

For the 'ones trick' to work you generally have to divide your new pmf by some large constant to ensure that the probability is between 0 and 1. While it appears that the pmf of the discrete weibull already ensures this, we have still added some large constant in the model anyways. So, here is what the model looks like now:

j.model ="
data{ 
C <- 10000
for(i in 1:N){
ones[i] <- 1
}
}
model{
for (i in 1:N){
discrete_weib[i] <- exp(-(y[i]/scale)^shape) - exp(-((y[i]+1)/scale)^shape)

ones[i] ~ dbern(discrete_weib[i]/C)
}

#priors
scale ~ dnorm(0, .0001) I(0, )
shape ~ dnorm(0, .0001) I(0, )
}
"

Note that we added 1) a vector of ones and a large constant in the data argument, 2) the pmf of the discrete weibull, and 3) we run that probability through a Bernoulli trial.

You can fit the model with the same code you have above, here is the summary which shows that the model successfully recovered the parameter values (scale = 7 and shape = 4).

       Lower95   Median  Upper95     Mean        SD Mode       MCerr MC%ofSD SSeff
scale 6.968277 7.289216 7.629413 7.290810 0.1695400   NA 0.001364831     0.8 15431
shape 3.843055 4.599420 5.357713 4.611583 0.3842862   NA 0.003124576     0.8 15126