# Negative binomial distributed variable with Poisson prior

I would like to know what is the unconditional distribution of variable $y \sim \text{Neg-Bin}(r, p)$ when $r \sim \text{Poisson}(\lambda)$.

• Are you sure you mean unconditional distribution..? – Tim Jan 30 '17 at 9:21
• Not sure about the term unconditional'. Sorry about that. I would need to know the distribution of $y$ given the Poisson prior on $r$. Thanks. – Simon74 Jan 30 '17 at 9:30
• @Tim My interpretation of the original question was that we have the conditional distribution of $Y|r$ which is $\text{NB}(r,p)$ and the distribution of $r$ which is $\text{Poisson}(\lambda)$. I think he wanted to find the distribution of $Y$ (or, as he says, the unconditional distribution), which would be $f_{Y}(y)=\sum_{x}f_{Y|X}(y|x)f_{X}(x)$. Not sure if that makes sense. I'm not sure the original question had a Bayesian context. – epp Jan 30 '17 at 11:08
• @StatsPlease this interpretation seems plausible if stated like this. Let's wait for OP's comment. Notice that he mentions prior. If this was the question, I'd edit my answer since it also applies to problem defined as so. Btw, I rolledback my edit. – Tim Jan 30 '17 at 11:13
• In a Bayesian context, if $y$ is data and $r$ a parameter, the 'unconditional', i.e., marginal distribution of $y$ would be interpreted as the prior predictive (which is what I also interpreted the question to be asking) – Juho Kokkala Jan 31 '17 at 5:58

As far as I know Poisson distribution is not a conjugate prior for $r$ parameter of negative binomial distribution. Nonetheless, if you need to estimate $r$ using Poisson prior, then it can be easily achieved using maximum a posteriori estimation using a direct search over a grid (i.e. you you estimate unnormalized posterior for all values on a grid of points and then find the maximum). Direct search is feasible in here since while in theory values of $r$ are countably infinite, in practice there would be a finite number of values with non-negligible probabilities, so you need to possibly check against $N+1$ values (integers $0,1,\dots,N$). If $r$ itself is not a huge value, then $N$ also does not need to be huge, so going through all the possible values would be very fast for any modern computer.

set.seed(123)
ptrue <- 0.32
rtrue <- 42
X <- rnbinom(200, rtrue, ptrue)

nmax <- 500
lambda <- 60

res <- numeric(nmax+1)
for (j in 0:nmax)
res[j+1] <- sum(dnbinom(X, j, ptrue, log = TRUE)) + dpois(j, lambda, log = TRUE)


As you can see on the plot below, using hyperparameter $\lambda=60$ it finds that the posterior distribution of $r$ reaches it's maximum at value $44$, that is pretty close to the true value equal to $42$.

plot(0:nmax, res, type = "l")
abline(v = which.max(res)+1, col = "red")
`

If you need to find the normalizing constant, then you also can evaluate the unnormalized posterior over a grid with some arbitrary large upper boundary (Poisson distribution has $[0, \infty)$ support) and sum them up, to obtain approximation of the Bayes theorem

$$g(r\mid x,p) = \frac{ \overbrace{f(x\mid r,p)}^\text{likelihood} \, \overbrace{g(r\mid\lambda)}^\text{prior} }{ \sum_{j=0}^{\infty} f(x\mid j,p) \, g(j\mid\lambda) }$$

where you would sum up to some arbitrary large integer instead of summing over all the possible values of $j$.

Honestly, even if analytical solution for this problem existed, then using a direct search (or other similar approaches) is so simple that it would less time to obtain estimates like this, then to find and verify the analytical solution.

Notice that the method described above could also be used if the question is understood like in the answer provided by StatsPlease. In such case for some $x$ of interest you would be summing over values of $r$ from zero up to some arbitrary large integer.

I think there is a little confusion around what exactly the question is asking. I know @Tim has provided an answer to this question, but my interpretation of the question was different.

As far as I can tell, you have the following distributions:

$$Y\,|\,r\sim \text{NB}(r,p)$$

and

$$r\sim \text{Poisson}(\lambda)$$

The distribution above, $Y\,|\,r$ is a conditional distribution. From what I can tell, you want to find the distribution of $Y$ not conditional on $r$. This isn't (as far as I'm aware) referred to as an unconditional distribution, but rather just the marginal distribution of $Y$.

Now, these distributions are discrete, and the expression for the density of $Y$ can be written as:

$$f_{Y}(y)=\sum_{r}f_{r,Y}(r,y)$$ where $f_{r,Y}(r,y)$ is the joint density of $r$ and $Y$. Now, this can be expressed as:

$$f_{Y}(y)=\sum_{r}f_{Y|r}(y|r)f_{r}(r)$$

Given that the densities are:

\begin{align} f_{Y|r}(y|r)&={{y+r-1}\choose{y}}(1-p)^{r}p^{y}\\ f_{r}(r)&=\frac{\lambda^{r}e^{-\lambda}}{r!} \end{align}

Now, I couldn't find a closed-form solution for the density of $Y$ (maybe someone else can). But I did run some simulations to get an idea for the distribution of $Y$ for a set of parameters $(p,\lambda)$. From the simulations run, I then found the MLE for the resulting data for a $\text{NB}(r,p)$ where $p$ is fixed (as it is a known parameter). The results were somewhat interesting:

$$Y\,|\,r\sim\text{NB}(r,p)\overset{*}{\Leftrightarrow} Y\sim\text{NB}(E[r]=\lambda,p)$$

$^{*}$This is purely speculative.

Essentially, the resulting $r$ parameter of the negative binomial distribution was always very close to:

$$r\approx E[r]=\lambda$$

Here is the comparison of the distribution obtained using the simulation results and the distribution where the negative binomial parameter $r=\lambda$. Once again, this is speculative, just thought it was interesting!

I looked a little deeper at this, for different values of $\lambda$, $p$ and simulations $N$.

For $\lambda=2$:

For $\lambda=10$:

For $\lambda=50$: