Vectorize Zero-truncated Poisson Density Computation

[Ben Lee]

Hello, 

In the past, I was able to speed up computation by vectorizing operations such as the density calculation for iid Bernoulli random variables:

dBernoulliVector <- nimbleFunction(

  run = function(x    = double(1),

                 prob = double(1),

                 log  = integer(0, default = 0)) { ## default should be zero

    returnType(double(0))

    logProb <- sum(dbinom(x, size = 1, prob = prob, log = TRUE))

    if(log) return(logProb) else return(exp(logProb))

  }

)

Is there a way for me to vectorize operations for Zero-truncated Poisson random variables?  I tried doing the following, but it doesn't seem to work:

dtpois <- nimbleFunction(

  run = function(x = integer(), lambda = double(), 

                 log = logical(0, default = 0)) {

    returnType(double())

    ## First handle non-zero data

    if(x==0){

      if(log){

        return(-Inf)

      }else{

        return(0)  

      }

      

    }else{

      if (log){

        return(dpois(x, lambda, log = TRUE) - log(1-exp(-lambda)))

      }else{

        return(dpois(x, lambda, log = FALSE)/(1-exp(-lambda)))

      }

    }

    

  })

rtpois <- nimbleFunction(

  run = function(n = integer(), lambda = double()) {

    returnType(integer())

    return(rpois(1, lambda)+1) # This was a lazy fix. I just added this component for completeness. I do not intend to sample from this distribution. 

  })

registerDistributions(list(

  dtpois = list(

    BUGSdist = "dtpois(lambda)",

    discrete = TRUE,

    range = c(0, Inf),

    types = c('value = integer()', 'lambda = double()')

  )))

dTruncPoissonVector <- nimbleFunction(

  run = function(x    = double(1),

                 lambda = double(1),

                 log  = integer(0, default = 0)) { ## default should be zero

    

    returnType(double(0))

    logProb <- sum(dtpois(x, lambda = lambda, log = TRUE))

    if(log) return(logProb) else return(exp(logProb))

  }

)

Thank you so much!

[Chris Paciorek]

Hi Ben,

It looks like your dtpois takes in a scalar but in dTruncPoissonVector you are passing in a vector to dtpois.

You can change dTruncPoissonVector to have a for loop that calls dtpois, if you want to keep the modularity you have, but my suggestion would be to have a single vectorized distribution that takes in vector x and does the truncated Poisson calculation in a vectorized fashion internally without calling out to a scalar distribution function. (That said, since this will be compiled, using a for loop within dTruncPoissVector that directly does the calculations you have in dtpois should be just as fast.)

-chris

--
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 on the web visit https://groups.google.com/d/msgid/nimble-users/111ab419-b00a-4535-b586-ec9950c8507dn%40googlegroups.com.