Modeling problem when N obeys Poisson distribution in binomial distribution

Dear All,

I now encounter a modeling problem, similar to assuming that the data obeys the binomial distribution, where N obeys the Poisson distribution.
y \sim binomial(N,p),\ \ \ \ N \sim poisson(lambda)

Since stan seems unable to handle the situation where N is positive infinity, I assume that there is an upper bound on N, such as 50.

I tried to write the model code

data {
  int<lower=0> N;
  int<lower=0> y[N];
  int<lower=0> max_N;

parameters {
  real<lower=0> lambda;
  real<lower=0, upper=1> p;
transformed parameters {
  vector[N] lp;
  for (i in 1:N){
    for  (j in (y[i]+1):max_N){
      lp[i] += poisson_lpmf(j | lambda) + binomial_lpmf(y[i] | j,p);

model {
  target += log_sum_exp(lp);

But got the following error.

Chain 4 Rejecting initial value:
Chain 4 Log probability evaluates to log(0), i.e. negative infinity.
Chain 4 Stan can’t start sampling from this initial value.

It seems that the calculation of which step got 0.

Thanks for any advice and it would help me a lot.


Unless, I am misunderstanding something about your model, it should be mathematically equivalent to just y \sim Poisson (p \lambda), so I would expect this form to greatly improve the behavior of your model.

More generally, to debug initial value rejected problems it is often helpful to use print statements to inspect intermediate values and see where the negative infinity creeps in.

Best of luck with your model!


Thank you for your reply

In fact, the model I really want to estimate is more complex.

\begin{equation} \begin{aligned} N &\sim possion(\lambda) \\ p &\sim zipf(s,N) = \frac{(1/k^s)}{\sum_{i=1}^{N}1/i^s}\\ q &= 1-(1-p)^M\\ n &\sim poisson\_binomial(q) \end{aligned} \end{equation}

Given M,n, estimation of \lambda and S

It is indeed a good point to transform the parameters.

I will try it again. Thank you very much!

What is k in your model ? The way you wrote this, it seems p is an integer-valued random variable, with support in [0, N].

Oops, sorry, I made a mistake in my expression.
p is a vector of length N. The value taken is determined by the zipf distribution.

So the second and third lines should probably be written as

\begin{equation} \begin{aligned} p_k&\sim zipf(k|s,N)\ \ (k=1,2,...,N)\\ q_k&=1-(1-p_k)^M \end{aligned} \end{equation}

Thank you for the correction!

sorry, I probably still don’t understand your model:

Is the Zipf distribution you refer to the same as what the Wikipedia has at Zipf's law - Wikipedia ? And does your Poisson Binomial match the formula at Poisson binomial distribution - Wikipedia ? If yes then p_k \geq 0 is an integer, right? But you need 0 < p_k < 1 for the Poisson binomial. But that is a contradiction, as with p_k with non-negative integer support, q_k can only be exactly 0, exactly 1 or smaller than zero… So I must be missing something.

Could you clarify where I am wrong?

1 Like