This is motivated by a problem where the count model is more complicated than the Poisson. But the general question is better understood with this simple model. Suppose I have a textbook hierarchical model:

X|Y \sim \text{Binomial}(Y, p),

Y \sim \text{Poisson}(\lambda).

Then it’s well known that X \sim \text{Poisson}(\lambda p). If you try and fit this model in Stan you’ll get a lot of divergences and the wrong parameter estimates, because only the product \lambda p is identified. The same problem crops up if you replace the Poisson with a negative binomial (under a mean/ inverse overdispersion parametrisation, `neg_binomial_2`

).

My question is: is there a way to reparametrise the model so as to identify \lambda and p?

Apologies in advance if this is a silly question.