# Sum of squared Negative Binomial probability masses

Let $(p_k)_{k=0, \dots, \infty}$ denote the probability masses of a Negative Binomial distribution with parameters $r>0$ and $p\in]0,1[$. I'm looking for the sum of their squares, $$\sum_{k=0}^\infty p_k^2$$ as a function of $r$ and $p$ (other parameterizations are also fine). In other words I am interested in (the exponential of) the second-order Renyi entropy of a Negative Binomial distribution.

Background:

• As per @whuber's comment here, this question is related, but not a duplicate... I'd still be interested in an analytical solution. – Stephan Kolassa Jul 21 '15 at 11:44
• I have edited your post to add the name second-order Renyi entropy. I think this could be useful for web search engines. I hope this is fine for you. – Stéphane Laurent Jul 21 '15 at 16:13

Take two independent Poisson random variables $X$ and $X'$, with means $\lambda$ and $\lambda'$. The formula answering your question in the Poisson case is a particular case of the identity $$\Pr(X = X' \mid \lambda, \lambda') = \exp(-\lambda)\exp(-\lambda')I_0(2\sqrt{\lambda}{\sqrt{\lambda'}}).$$

> lambda <- 1; lambdaa <- 2
> sum(dpois(0:100,lambda)*dpois(0:100,lambdaa))
[1] 0.2117121
> gsl::bessel_I0(2*sqrt(lambda*lambdaa)) * exp(-lambda-lambdaa)
[1] 0.2117121


Your problem is the same as calculating $\int\int\Pr(X= X' \mid \lambda, \lambda') f_{a,b}(\lambda)f_{a,b}(\lambda')d\lambda d\lambda'$ where $f_{a,b}$ is the $\Gamma(a,b)$ pdf, setting $r=a$ and $p=\frac{b}{1+b}$ with your notations, because of the link between the negative binomial distribution and the Poisson-Gamma distribution.

Let's start by $$\int \exp(-\lambda)I_0(2\sqrt{\lambda}{\sqrt{\lambda'}})f_{a,b}(\lambda)d\lambda = \frac{b^a}{\Gamma(a)}\int \lambda^{a-1}\exp\bigl(-(b+1)\lambda\bigr)I_0(2\sqrt{\lambda}{\sqrt{\lambda'}})d\lambda.$$

According to Mathematica this is equal to $${\left(\frac{b}{1+b}\right)}^a {}_1\!F_1\left(a, 1, \frac{\lambda'}{b+1}\right)$$ where ${}_1\!F_1$ is the Kummer hypergeometric function.

Now we can even get something for $$\begin{multline} \int \exp(-\lambda'){}_1\!F_1\left(a, 1, \frac{\lambda'}{b+1}\right)f_{a',b'}(\lambda')d\lambda' \\ = \frac{{b'}^{a'}}{\Gamma(a')}\int {\lambda'}^{a'-1}\exp\bigl(-(b'+1)\lambda'\bigr) {{}_1\!F_1}\left(a, 1, \frac{\lambda'}{b+1}\right)d\lambda'. \end{multline}$$ Indeed, Mathematica gives $${\left(\frac{b'}{1+b'}\right)}^{a'} {}_2\!F_1\left(a, a', 1, \frac{1}{(b+1)(b'+1)}\right)$$ where ${}_2\!F_1$ is the Gauss hypergeometric function. The final result is beautiful: $${\left(\frac{b}{1+b}\right)}^{a}{\left(\frac{b'}{1+b'}\right)}^{a'}{}_2\!F_1\left(a, a', 1, \frac{1}{(b+1)(b'+1)}\right),$$ and even a bit more beautiful with your notations: $$p^{a}{p'}^{a'}{}_2\!F_1\left(a, a', 1, (1-p)(1-p')\right)$$

Check:

> a <- 2; A <- 3; b <- 5; B <- 8
> (b/(1+b))^a*(B/(1+B))^A*gsl::hyperg_2F1(a,A,1,1/(b+1)/(B+1))
[1] 0.5450618
> sum(dnbinom(0:100, a, b/(1+b))*dnbinom(0:100, A, B/(1+B)))
[1] 0.5450618


In the special case you are interested in, the sum of squares is $$p^{2a}{}_2\!F_1\left(a, a, 1, {(1-p)}^2\right),$$ and the second-order Renyi entropy is $$-2a \log p - \log {}_2\!F_1\left(a, a, 1, {(1-p)}^2\right).$$

• Thank you. This looks both intriguing and intimidating. I'll need a little time to digest it. In the meantime, could you please introduce me to your friend ${}_kF_\ell$? – Stephan Kolassa Jul 21 '15 at 15:28
• @StephanKolassa I have edited my answer to add their names. You will easily find some information about them. – Stéphane Laurent Jul 21 '15 at 15:54