I'm trying to implement the following model in Stan: $$\begin{align} \text{Pr}(y|n,p) & \sim \text{Binomial}(n,p)\\ \text{Pr}(n|\lambda) & \sim \text{Poisson}(\lambda) \end{align}$$
In this model, $y$ and $n$ are non-negative integers, $0<p<1$, $\lambda > 0$, and $y \le n$. Because Stan does not allow for discrete parameters, I need to marginalize over $n$ (i.e., find the probability mass function of $y$ given $p$ and $\lambda$. The PMFs of $y$ and $n$ are as follows:
$$\begin{align} \text{Pr}(y|n,p) = & \frac{n!}{y!(n-y)!}p^{y}(1-p)^{n-y}\\ \text{Pr}(n|\lambda) =& \frac{\lambda^{n}e^{-\lambda}}{n!}\end{align}$$
My attempt to marginalize over $n$ is as follows:
$$\begin{alignat}{4} \text{Pr}(y|p,\lambda) & = \sum_{n=y}^{\infty}\text{Pr}(y|n,p)\text{Pr}(n|\lambda)\\\\ & = \sum_{n=y}^{\infty}\frac{\lambda^{n}e^{-\lambda}}{y!(n-y)!}p^{y}(1-p)^{n-y}\\\\ & = \frac{e^{-\lambda p}p^{y}\lambda^{y}}{\Gamma(y+1)} \end{alignat} $$
The final step was done by Wolfram Alpha, as my own mathematical skills are a bit rusty. Does this seem correct?
To implement this, I first took the log PMF:
$$\text{log-Pr}(y|p,\lambda)=-\lambda p+y(\text{log}p+\text{log}\lambda)-\text{log}\Gamma(y+1)$$
The Stan code for the function is:
real binomial_poisson_lpmf(int[] y, vector p, vector lambda) { // vectorized; assumes that y, p, and lambda are all of the same length.
vector[size(y)] out;
for(i in 1:size(y))
out[i] = -lambda[i] * p[i] + y[i] * (log(p[i]) + log(lambda[i])) - lgamma(y[i]+1));
return(sum(out));
}
Anyway, I would appreciate any feedback on whether or not this seems valid; I would also like to know if anyone knows of any citations on the subject.
Thank you for your time.