Vectorized Poisson Error

[Colton Padilla]

Hi everyone,

This is a relatively simple question and hopefully is a simple answer. I am trying to sample using a vectorized poisson to see if it will speed my model up. I have gotten the distribution function correct but it seems that I have the random generation function incorrect. Could someone look this code over and see if there are any errors?

# Vectorized poisson
dpois_vec <- nimbleFunction(
  run = function(x = double(1),
                 lambda = double(1),
                 log = integer(0, default = 0)){
    returnType(double(0))
    logProb <- sum(dpois(x, lambda, log = TRUE))
    if(log) return(logProb)
    else return(exp(logProb))
  }
)

# Random generation
rpois_vec <- nimbleFunction(
  run = function(n = integer(0, default = 1),
                 lambda = double(1)){
    returnType(double(1))
    return(rpois(length(lambda), lambda))
  }
)

[Chris Paciorek]

I don't see anything wrong. Can you show us what the error is that you are getting?

--
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/e241be44-cdcb-475f-8810-ee88a9a35bf1n%40googlegroups.com.

[Colton Padilla]

It is telling me that the logProb of the data node is NA or NaN. Here is the code of what I am trying to do. The dec_pop variable is the output of a population projection from an ODE. I am able to easily run the model with the non-vectorized version shown at the bottom of the code, but when I try to vectorize, I keep getting the value/logProb is NA or NaN. 

# Vectorized process model
  J_est[1:maxtime] ~ dpois_vec(dec_pop[1:maxtime, 1])
  Y_est[1:maxtime] ~ dpois_vec(dec_pop[1:maxtime, 2])
  A_est[1:maxtime] ~ dpois_vec(dec_pop[1:maxtime, 3])
 
  # Vectorized observation model
  J[2:maxtime] ~ dpois_vec(J_est[2:maxtime])
  Y[2:maxtime] ~ dpois_vec(Y_est[2:maxtime])
  A[2:maxtime] ~ dpois_vec(A_est[2:maxtime])
 
  # Non-vectorized process model
  for(t in 1:maxtime){
  J_est[t] ~ dpois(dec_pop[t, 1])
  Y_est[t] ~ dpois(dec_pop[t, 2])
  A_est[t] ~ dpois(dec_pop[t, 3])
  
  # Observation model
  for(t in 2:maxtime){
   J[t] ~ dpois(J_est[t])
   Y[t] ~ dpois(Y_est[t])
   A[t] ~ dpois(A_est[t])
  

[Colton Padilla]

I can also do the vectorized process model and then loop the observation model, so I am guessing it has something to do with the random number generation function.

[Chris Paciorek]

If you just compile `rpois_vec` and run it, it should work fine - I just tried it.

So I think the question is why you're getting unexpected behavior in the model. It may have to do with model initialization.

I would try using `simulate` directly with the compiled model and see what happens in terms of the values of J, Y, A, etc. and their logProb values. On a related note, we have suggestions here on how to debug model initialization.

-chris

To view this discussion on the web visit https://groups.google.com/d/msgid/nimble-users/4714fbb7-fda7-4eec-916b-dea6d55e83b4n%40googlegroups.com.