Random Poisson variable is a negative but must be nonnegative!

Modeling
,
1

I am getting started with Stan, and I have the following state-process model coded in Stan:

data {
  int<lower=0> maxIter;
  int<lower=0> numGears;
  int C[numGears,maxIter,2];
  real totalArea;
  real a[numGears];
}

parameters {
  real<lower=0> N0;
  real<lower=0,upper=1> juvProp;
  real<lower=0> f;
  real<lower=0,upper=1> sJ;
  real<lower=0,upper=1> b;
  real<lower=0,upper=1> sA;
  real<lower=0,upper=1> c[numGears];
}

transformed parameters{
  real J[maxIter];
  real A[maxIter];
  real B[maxIter];
  
  J[1] = N0*juvProp;
  A[1] = N0*(1-b)*(1-juvProp);
  B[1] = N0*b*(1-juvProp);
  for(i in 1:(maxIter-1)){
    J[i+1] = f*B[i]/2;
    A[i+1] = (1-b)*sJ*J[i]+sA*A[i];
    B[i+1] = b*sJ*J[i];
  }
}

model {
  int nA[numGears,maxIter];
  int nJ[numGears,maxIter];
  
  for(i in 1:numGears){
    c[i] ~ uniform(0,1);
		for(j in 1:maxIter){
		  nJ[i,j] ~ poisson(a[i]*J[j]/totalArea);
			C[i,j,1] ~ binomial(nJ[i,j],c[i]);
			nA[i,j] ~ poisson(a[i]*(A[j]+B[j])/totalArea);
			C[i,j,2] ~ binomial(nA[i,j],1-c[i]);
		}
	}
	
	N0 ~ normal(1000,10000) T[0,];
	juvProp ~ uniform(0,1);
	f ~ lognormal(0,5);
	sJ ~ uniform(0,1);
	b ~ uniform(0,1);
	sA ~ uniform(0,1);
}

When I try to fit the model in rstan:

library(rstan)
numGears <- 2
maxIter <- 10
a <- c(0.2,0.3)
totalArea <- 1
C <- array(rbinom(1000,size = 1000,0.3),c(2,10,2))
stanData <- list(maxIter=maxIter,numGears=numGears,C=C,totalArea=totalArea,a=a)

stanModel <- stan(
  file = "StanModel.stan",
  data = stanData,
  chains= 1,
  warmup = 5000,
  iter = 50000,
  cores = 8
)

I get the following error:

SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 1).
Chain 1: Rejecting initial value:
Chain 1:   Error evaluating the log probability at the initial value.
Chain 1: Exception: poisson_lpmf: Random variable is -2147483648, but must be nonnegative! (in 'string', line 38, column 4 to column 43)

This doesn’t make much sense to me. All the parameters are set to be positive, so nothing should produce a negative number during sampling.

What am I doing wrong?

2

You are expecting sampling statements like nJ[i,j] ~ poisson(a[i]*J[j]/totalArea); to populate nJ[i,j] with a value sampled from a Poisson distribution. But that’s not what these statements achieve in Stan. Instead, in Stan, these statements exist to define something proportional to the joint log-density (or log-mass) of the data and the parameters. I’m not too clear on why nJ[i,j] is getting initialized as -2147483648, but what you really need to do for latent parameters is to declare them as parameters, so that they will get initialized to some appropriate value subject to the appropriate constraints, and proceed from there.

However, this brings up a separate problem, because you cannot have discrete parameters in Stan. So to fit this model, what you really need is to marginalize nA and nJ out of the likelihood. Fortunately, these poisson-binomials marginalize straightforwardly to poissons, so you can just write, for example C[i,j,1] ~ poisson(c[i] * a[i] * J[j] / totalArea)

Note that this marginalization works only because each sampled poisson variate is only used once in the binomial likelihood. In cases where a Poisson sample serves as the number of counts in multiple binomial samples (e.g. in N-mixture models), you need an alternative scheme marginalize over the variate, e.g. by brute-force summing through the terms up to some bound of vanishingly low upper-tail probability in the Poisson.

2 Likes
3

A post was split to a new topic: N-mixture likelihoods

4

Thanks that is a great point! If it were the case that the math didn’t work out as in this case, would an alternative solution be to include nJ and nA as dummy variables in the data section and let Stan deal with them? Or are data mutations not allowed?

5

Data mutations are allowed in transformed data, but these must be deterministic data transformations that run once prior to sampling. You can never use latent discrete parameters in your likelihood computation in Stan; these must always be marginalized out.

1 Like
6

Thanks! I’m guessing that in that case one would need to do something like what you mentioned in the other post.

7

yes, that’s right :)