Build distribution for generalized Poisson

[Xinying Fang]

Hi,

I'd like to build an under-dispersion Poisson model with generalized Poisson distribution. The PMF formula is attached in the email. Here are my code to build the distribution in NIMBLE:

## ---- GP-1 log PMF (density) ----
dgenpois_gp1 <- nimbleFunction(
  run = function(x   = integer(0),
                 mu  = double(0),
                 phi = double(0),
                 log = integer(0, default = 0)) {
    returnType(double(0))          # <-- REQUIRED
   
    # constraints
    if (mu <= 0 | phi >= 1 | x < 0) {
      if (log) return(-Inf) else return(0.0)
    }
   
    lambda <- mu * (1 - phi)
    base   <- lambda + phi * x
   
    if (base <= 0) {
      if (log) return(-Inf) else return(0.0)
    }
   
    logProb <- log(lambda) + (x - 1) * log(base) - base - lgamma(x + 1.0)
    if (log) return(logProb) else return(exp(logProb))
  }
)

## ---- GP-1 RNG ----
rgenpois_gp1 <- nimbleFunction(
  run = function(n   = integer(0),
                 mu  = double(0),
                 phi = double(0)) {
    returnType(integer(0))         # <-- REQUIRED
   
    if (n != 1 | mu <= 0 | phi >= 1) return(NA_integer_)
   
    lambda <- mu * (1 - phi)
    u <- runif(1, 0, 1)
    y <- 0L
    cdf <- 0.0
   
    # finite support when phi < 0
    yMax <- if (phi < 0) floor(lambda / (-phi) - 1e-12) else 10000000L
   
    while (y <= yMax) {
      base <- lambda + phi * y
      if (base <= 0) {
        # exit the loop without 'break'
        y <- yMax + 1L
      } else {
        # accumulate probability
        # log pmf: log(lambda) + (y-1)*log(base) - base - lgamma(y+1)
        cdf <- cdf + exp(log(lambda) + (y - 1) * log(base) - base - lgamma(y + 1.0))
        if (u <= cdf) return(y)
        y <- y + 1L
      }
    }
    return(as.integer(yMax))
  }
)

## ---- Register the distribution with NIMBLE ----
registerDistributions(list(
  dgenpois_gp1 = list(
    BUGSdist = "dgenpois_gp1(mu, phi)",
    types    = c("value = integer()", "mu = double()", "phi = double()"),
    pqAvail  = FALSE
  )
))

When I perform a test on the model:

set.seed(1)
N  <- 200
x  <- rnorm(N)
beta0 <- 0.2
beta1 <- 0.8
phi_true <- -0.3                 # underdispersed example (Var < mean)
mu  <- exp(beta0 + beta1 * x)

gp1_code <- nimbleCode({
  for (i in 1:N) {
    y[i] ~ dgenpois_gp1(mu[i], phi)
    log(mu[i]) <- beta0 + beta1 * x[i]
  }
  # diffuse priors
  beta0 ~ dnorm(0, 1/100)     # precision parameterization
  beta1 ~ dnorm(0, 1/100)
  # keep phi safely away from 1; allow under/over-dispersion
  phi ~ dunif(-0.9, 0.95)
})

constants <- list(N = N, x = x)
data      <- list(y = y)
inits     <- list(beta0 = 0, beta1 = 0, phi = 0)

gp1_model <- nimbleModel(code = gp1_code, constants = constants,
                         data = data, inits = inits)

Cmodel <- compileNimble(gp1_model)

I met this error:

> Cmodel <- compileNimble(gp1_model)   
Compiling
  [Note] This may take a minute.
  [Note] Use 'showCompilerOutput = TRUE' to see C++ compilation details.
Error in code$args[[i]]$isCall : $ operator is invalid for atomic vectors

Could you help me check is there anything wrong with the built distribution?

Thank you so much!

Best,

Xinying

[Perry de Valpine]

Hi Xinying,

Thanks for the question.

You can try compiling the nimbleFunctions separately to isolate the problem.

Cdgenpois_gp1 <- compileNimble(dgenpois_gp1) # compiles without error. I didn't test if it is correct.

Crgenpois_gp1 <- compileNimble(rgenpois_gp1) # gives the error you saw when compiling the model that uses this distribution.

I see in rgenpois_gp1, there are several things that would work in R but are not supported for compilation in nimbleFunctions.

Model variables are always doubles, so the returnType should be double(), and neither "as.integer" nor NA_integer_ are supported.

Functions supported for compilation are in the User Manual at r-nimble.org (which is up again).

if-then-else does not return a value in compiled nimbleFunctions.

The following version compiles. Please see if it works in your model.

rgenpois_gp1 <- nimbleFunction(
  run = function(n   = integer(0),
                 mu  = double(0),
                 phi = double(0)) {

    returnType(double(0))         # <-- REQUIRED
   
    if (n != 1 | mu <= 0 | phi >= 1) return(NA)

   
    lambda <- mu * (1 - phi)
    u <- runif(1, 0, 1)
    y <- 0L
    cdf <- 0.0
   
    # finite support when phi < 0

if(phi < 0)

  yMax <- floor(lambda/(-phi) -1e-12)

else

  yMax <- 10000000 # yMax can't be a double in one place and an integer somewhere else. Types are static as in C++.
#    yMax <- if (phi < 0) floor(lambda / (-phi) - 1e-12) else 10000000L

   
    while (y <= yMax) {
      base <- lambda + phi * y
      if (base <= 0) {
        # exit the loop without 'break'
        y <- yMax + 1L
      } else {
        # accumulate probability
        # log pmf: log(lambda) + (y-1)*log(base) - base - lgamma(y+1)
        cdf <- cdf + exp(log(lambda) + (y - 1) * log(base) - base - lgamma(y + 1.0))
        if (u <= cdf) return(y)
        y <- y + 1L
      }
    }

    return(yMax)
  }
)

In registerDistributions, the value should be double() (and the automatic registration might work fine).

HTH

Perry

--
You received this message because you are subscribed to the Google Groups "nimble-users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to nimble-users...@googlegroups.com.
To view this discussion visit https://groups.google.com/d/msgid/nimble-users/8c031601-2357-46d7-a5bf-ede4e39d5fb9n%40googlegroups.com.