# Sum of Binomial and Poisson random variables

If we have two independent random variables $X_1 \sim \mathrm{Binom}(n,p)$ and $X_2 \sim \mathrm{Pois}(\lambda)$, what is the probability mass function of $X_1 + X_2$?

NB This is not homework for me.

• I guess you tried convoluting? en.wikipedia.org/wiki/… Where did you get stuck? I assume there is no closed form, otherwise the solution would probably be here: en.wikipedia.org/wiki/… Nov 28, 2012 at 12:03
• Yes that's what I tried, but maybe I have found an answer here: mathstatica.com/SumBinomialPoisson Kummer confluent hypergeometric function..hugh! Nov 28, 2012 at 12:07
• I've readded the homework tag in accordance with its use on this site. Cheers. :-) Nov 28, 2012 at 13:35
• Novel means new (not known or published before). I also don't agree that using known methods to solve new problems makes it homework -- the same can be said for the majority of journal articles publishing results on distributions. Apr 22, 2019 at 5:06
• As in many other cases in statistics where a hypergeometric function appears with integral arguments, you can understand it to be a shorthand notation for the implicit (finite) sum in the convolution if you wish. The advantage of such an expression is that there are myriad ways to manipulate it into simpler forms and it can often be evaluated without actually performing the summation.
– whuber
Apr 23, 2019 at 2:22

You will end up with two different formulas for $p_{X_1+X_2}(k)$, one for $0 \leq k < n$, and one for $k \geq n$. The easiest way of doing this problem is to compute the product of $\sum_{i=0}^n p_{X_1}(i)z^k$ and $\sum_{j=0}^{\infty}p_{X_2}(j)z^j$. Then, $p_{X_1+X_2}(k)$ is the coefficient of $z^k$ in the product. No simplification of the sums is possible.

Giving the closed formula in terms of generalized hypergeometric functions (GHF) hinted at in other answers (the GHF in this case is really only a finite polynomial, so is a shorthand for the finite sum.) I used maple to sum the convolution, with this result: $$\DeclareMathOperator{\P}{\mathbb{P}} \P(X_1+X_2=k)= \sum_{x_1=0}^{\min(n,k)} \binom{n}{x_1} p^{x_1}(1-p)^{n-x_1} e^{-\lambda} \frac{\lambda^{k-x_1}}{(k-x_1)!}= {\frac { \left( 1-p \right) ^{n}{{\rm e}^{-\lambda}}{\lambda}^{k}}{ \Gamma \left( k+1 \right) } {\mbox{_2F_0}(-k,-n;\,\ ;\,-{\frac {p}{ \left( p-1 \right) \lambda}})} }$$

Dilip Sarwate stated 7 years ago that no simplification is possible, although this has been challenged in comments. However, I think it is useful to note that even without any simplification the computation is quite straightforward in any spreadsheet or programming language.

Here is an implementation in R:

# example parameters
n <- 10
p <- .3
lambda <- 5

# probability for just a single value
x <- 10  # example value
sum(dbinom(0:x, n, p) * dpois(x:0, lambda))

# probability function for all values
x0  <- 0:30   # 0 to the maximum value of interest
x   <- outer(x0, x0, "+")
db  <- dbinom(x0, n, p)
dp  <- dpois(x0, lambda)
dbp <- outer(db, dp)
aggregate(as.vector(dbp), by=list(as.vector(x)), sum)[1:(max(x0)+1),]

• Dilip did not show that no simplification of the sums is possible: he stated such an assertion (and the assertion does not appear to be correct). If you follow the links provided by the OP, a solution is provided in terms of Kummer confluent hypergeometric functions. Apr 21, 2019 at 16:17
• @wolfies - That would be a very interesting point in a new answer to this old question. Probably more interesting than mine.
– Pere
Apr 21, 2019 at 16:21
• A potentially faster approach for large n in the binomial, and large lambda would involve fast Fourier transforms (or similar). I've successfully used it on a number of real-world problems where the convolution is not algebraically convenient, but numerical answers are sufficient, and where multiple independent variates were being added. Apr 24, 2019 at 6:03
• Re @Glen_b's comment, for larger values of $n$ and $\lambda$ this brute-force convolution becomes cumbersome. Moreover, the challenge is not to implement it, but to find suitable endpoints for computing the dpois array: fixing x at 10 obviously won't cut it. One reliable method is to set x to extreme percentiles of the distribution, such as x<-qpois(0:1+c(1,-1)*1e-6, lambda), then compute dpois for the range of x, and then "chop" the results (with zapsmall) before proceeding with the outer product. When n is large, apply a similar procedure to the binomial probabilities.
– whuber
Apr 24, 2019 at 14:30
• Indeed. I did something similar with my own application -- going out sufficiently far gave the required quantiles as accurately as were needed. Apr 25, 2019 at 4:25