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.
Cross Validated is a question and answer site for people interested in statistics, machine learning, data analysis, data mining, and data visualization. It only takes a minute to sign up.
Sign up to join this communityIf 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.
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{$_2$F$_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),]
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.
$\endgroup$
– whuber♦
Apr 24 '19 at 14:30