Using wikipedia I found a way to calculate the probability mass function resulting from the sum of two Poisson random variables. However, I think that the approach I have is wrong.

Let $X_1, X_2$ be two independent Poisson random variables with mean $\lambda_1, \lambda_2$, and $S_2 = a_1 X_1+a_2 X_2$, where the $a_1$ and $a_2$ are constants, then the probability-generating function of $S_2$ is given by $$ G_{S_2}(z) = \operatorname{E}(z^{S_2})= \operatorname{E}(z^{a_1 X_1+a_2 X_2}) G_{X_1}(z^{a_1})G_{X_2}(z^{a_2}). $$ Now, using the fact that the probability-generating function for a Poisson random variable is $G_{X_i}(z) = \textrm{e}^{\lambda_i(z - 1)}$, we can write the probability-generating function of the sum of the two independent Poisson random variables as $$ \begin{aligned} G_{S_2}(z) &= \textrm{e}^{\lambda_1(z^{a_1} - 1)}\textrm{e}^{\lambda_2(z^{a_2} - 1)} \\ &= \textrm{e}^{\lambda_1(z^{a_1} - 1)+\lambda_2(z^{a_2} - 1)}. \end{aligned} $$ It seems that the probability mass function of $S_2$ is recovered by taking derivatives of $G_{S_2}(z)$ $\operatorname{Pr}(S_2 = k) = \frac{G_{S_2}^{(k)}(0)}{k!}$, where $G_{S_2}^{(k)} = \frac{d^k G_{S_2}(z)}{ d z^k}$.

Is this is correct? I have the feeling I cannot just take the derivative to obtain the probability mass function, because of the constants $a_1$ and $a_2$. Is this right? Is there an alternative approach?

If this is correct can I now obtain an approximation of the cumulative distribution by truncating the infinite sum over all k?

  • 1
    $\begingroup$ Some particular kinds of combinations have names (e.g. the Skellam distribution), but the general case its just compound Poisson, I think. There will be no closed form expression that's simpler than the original sum in general, though I suppose you can just write it as a convolution. When you say 'generate' do you mean 'simulate random values from', or something else? $\endgroup$
    – Glen_b
    Commented Feb 8, 2013 at 2:56
  • $\begingroup$ A closely related special case--that of the average of Poisson variates, binned by integers--is treated at stats.stackexchange.com/questions/35042/…. $\endgroup$
    – whuber
    Commented Feb 8, 2013 at 3:17
  • 1
    $\begingroup$ Why are you scaling the summands with $a_1$ and $a_2$? The sum is just another Poisson distribution without this. The variables take values in the positive integers, so something like $1$ times the first plus $\sqrt{2}$ times the second is usually quite unnatural, and would let you recover the values of both variables. $\endgroup$ Commented Feb 9, 2013 at 19:57
  • 4
    $\begingroup$ The difficulty here is that unless both $a_1$ and $a_2$ are integers, one cannot be sure that $S_2$ takes on integer values only. Thus, you need to find not just $P(S_2 = k)$ for integer values of $k$ but also $P(S_2 = \alpha)$ for each $\alpha$ that can be expressed as $a_1m + a_2n$ for nonnegative integers $m$ and $n$. $\endgroup$ Commented Feb 9, 2013 at 21:30
  • 1
    $\begingroup$ I don't think you can do much better than a brute-force approach which finds the possible values that $S_2$ can take on and then for each $\alpha$, use $$P\{S_2 = \alpha\} = \sum_{a_1m + a_2n = \alpha}P\{X_1=m\}P\{X_2=n\} = \sum_{a_1m + a_2n = \alpha} \exp(-\lambda_1m)\frac{\lambda_1^m}{m!}\exp(-\lambda_2n)\frac{\lambda_2^n}{n!}.$$ For most choices of $a_1$ and $a_2$, I would expect that most sums will reduce to a single term. I expect you know that for $a_1=a_2=1$, $S_2$ is a Poisson random variable with parameter $\lambda_1+\lambda_2$. $\endgroup$ Commented Feb 9, 2013 at 22:01

3 Answers 3


Provided not a whole lot of probability is concentrated on any single value in this linear combination, it looks like a Cornish-Fisher expansion may provide good approximations to the (inverse) CDF.

Recall that this expansion adjusts the inverse CDF of the standard Normal distribution using the first few cumulants of $S_2$. Its skewness $\beta_1$ is

$$\frac{a_1^3 \lambda_1 + a_2^3 \lambda_2}{\left(\sqrt{a_1^2 \lambda_1 + a_2^2 \lambda_2}\right)^3}$$

and its kurtosis $\beta_2$ is

$$\frac{a_1^4 \lambda_1 + 3a_1^4 \lambda_1^2 + a_2^4 \lambda_2 + 6 a_1^2 a_2^2 \lambda_1 \lambda_2 + 3 a_2^4 \lambda_2^2}{\left(a_1^2 \lambda_1 + a_2^2 \lambda_2\right)^2}.$$

To find the $\alpha$ percentile of the standardized version of $S_2$, compute

$$w_\alpha = z +\frac{1}{6} \beta _1 \left(z^2-1\right) +\frac{1}{24} \left(\beta _2-3\right) \left(z^2-3\right) z-\frac{1}{36} \beta _1^2 z \left(2 z^2-5 z\right)-\frac{1}{24} \left(\beta _2-3\right) \beta _1 \left(z^4-5 z^2+2\right)$$

where $z$ is the $\alpha$ percentile of the standard Normal distribution. The percentile of $S_2$ thereby is

$$a_1 \lambda_1 + a_2 \lambda_2 + w_\alpha \sqrt{a_1^2 \lambda_1 + a_2^2 \lambda_2}.$$

Numerical experiments suggest this is a good approximation once both $\lambda_1$ and $\lambda_2$ exceed $5$ or so. For example, consider the case $\lambda_1 = 5,$ $\lambda_2=5\pi/2,$ $a_1=\pi,$ and $a_2=-2$ (arranged to give a zero mean for convenience):


The blue shaded portion is the numerically computed CDF of $S_2$ while the solid red underneath is the Cornish-Fisher approximation. The approximation is essentially a smooth of the actual distribution, showing only small systematic departures.

  • 2
    $\begingroup$ Nice use of an often-forgotten tool... and, of course, for either of $\lambda_1$ or $\lambda_2 \leq 5$ or so, the brute force convolution method won't be all that painful. $\endgroup$
    – jbowman
    Commented Oct 17, 2013 at 16:30

Use the convolution:

Let $f_{X_1}(x_1)= \dfrac{\lambda^{x_1}e^{-\lambda}}{x_1!} $ for $x_1 \geq 0$, $f_{X_1}(x_1)= 0$ otherwise, and $f_{X_2}(x_2)=\dfrac{\lambda^{x_2}e^{-\lambda}}{x_2!} $ for $x_2 \geq 0$, $f_{X_2}(x_2)= 0$ otherwise.

Let $Z=X_1+X_2\rightarrow X_1=Z-X_2$, so $$f_Z(z)=\int\limits_{-\infty}^{\infty}\int\limits_{-\infty}^{\infty}f_{x_1,x_2}(z-x_2,x_2)dx_1dx_2$$ The former is known as convolution.

If $X_1$ and $X_2$ are independent, $$f_Z(z)=\int\limits_{-\infty}^{\infty}\int\limits_{-\infty}^{\infty}f_{X_1}(z-x_2)f_{X_2}(x_2)dx_1dx_2$$ This way you can obtain the distribution of the sum of two continuous random variables.

For the discrete poisson distribution $$f_Z(z)=\sum\limits_{x_2=0}^{z} \dfrac{\lambda^{z-x_{2}}_1e^{-\lambda_1}}{(z-x_2)!}\dfrac{\lambda^{x_2}_2e^{-\lambda_2}}{x_2!}$$ $$= e^{-(\lambda_1+\lambda_2)}\dfrac{(\lambda_1+\lambda_2)^z}{z!}$$ Which is also a Poisson distribution with parameter $\lambda_1+\lambda_2$

  • 3
    $\begingroup$ This appears to answer a different question: namely, how to add two Poisson distributions. It is the special case $a_1=a_2=1$ (but can be extended to the cases $a_1 = a_2$ without any trouble). But what would you do when $a_1 \ne a_2$? $\endgroup$
    – whuber
    Commented Oct 17, 2013 at 17:00

I think the solution is the concept of a compound Poisson distribution. The idea is a random sum $$ S = \sum_{i=1}^N X_i $$ with $N$ Poisson distributed and $X_i$ and $iid$ sequence independent of $N$. When we restric to the case that $X_i=k$ always, then we can describe $k N$ for a real number $k$ and a Poisson distributed $N$. You get the pgf by $$ E[s^{k N}] = E[(s^{k})^N] = G_N(s^{k}) = \exp(\lambda(s^k-1)) $$ For the sum $Z = k_1 N_1 + k_2 N_2$ you get $$ G_Z(s) = \exp(\lambda_1(s^{k_1}-1) + \lambda_2(s^{k_2}-1)). $$ define $\lambda = \lambda_1 + \lambda_2$ then $$ G_Z(s) = \exp(\lambda ( \frac{\lambda_1}{\lambda}(s^{k_1}-1)+ \frac{\lambda_2}{\lambda}(s^{k_1}-1)) = \exp(\lambda (\frac{\lambda_1}{\lambda}s^{k_1}+ \frac{\lambda_2}{\lambda}s^{k_1}-1)). $$ The final interpretation is that the resulting rv is a compound Poisson distribution with intensity $\lambda = \lambda_1 + \lambda_2$ and distribution of the $X_i$ that take the value $k_1$ with probability $\lambda_1/\lambda$ and the value $k_2$ with $\lambda_2/\lambda$.

Having proved that the distributions is compound Poisson we can either use Panjer recursion in the case that $k_1$ and $k_2$ are positive integers. Or we can easily derive the Fourier transform from the form of the pgf and get the distribution back by the inverse. Note that there is a point mass at $0$.

Edit after a discussion:

I think the best you can do is MC. You could use the derivation that this is a compound Poisson distr.

  1. sample N from $Pois(\lambda)$ (very efficient)
  2. then for each $i=1,\ldots,N$ sample whether it is from $X_1$ or $X_2$ where the probability of the first is $\lambda_1/\lambda$. Do this by sampling a Bernoulli rv with probability of success $\lambda_1/\lambda$. If it is $1$ then add $k_1$ to the sampled sum else add $k_2$.

You will have a sample of say 100 000 in seconds.

Alternatively you can sample the two summands in your inital representation separately ... this will be as quick.

Everything else (FFT) is complicated if the constant factor k1 and k2 are totally general.

  • 1
    $\begingroup$ And the final distribution can be found by the Panjer algorithm if the factors are integers. $\endgroup$
    – Richi W
    Commented Feb 10, 2013 at 0:13
  • $\begingroup$ Thanks! I got to $G_{S_2}(z) = \textrm{e}^{\lambda_1(z^{a_1} - 1)}\textrm{e}^{\lambda_2(z^{a_2} - 1)}$ However, starting from this I would like to find a way to be able to obtain some sort of distribution. You mentioned the Panjer algorithm? However, in this case $a_1,a_2 \in R^1$. @DilipSarwate Just mentioned that it is impossible to simplify the following $P\{S_2 = \alpha\} = \sum_{a_1m + a_2n = \alpha}P\{X_1=m\}P\{X_2=n\} = \sum_{a_1m + a_2n = \alpha} \exp(-\lambda_1m)\frac{\lambda_1^m}{m!}\exp(-\lambda_2n)\frac{\lambda_2^n}{n!},$ for generally $a_1,a_2 $. $\endgroup$
    – Michel
    Commented Feb 10, 2013 at 3:13
  • $\begingroup$ Hi Michel, I edited my response. Yes Panjer is of limited use. But you could try a Fourier transform approach. However non-integer units are problematic ... I have to think more about what to do in this case. Either way it is important to note that the result is a compound Poisson distribution (not a "simple" Poisson distribution). $\endgroup$
    – Richi W
    Commented Feb 10, 2013 at 13:19
  • $\begingroup$ Hi Richard, Thanks for your update! Do you mean that I should numerically try to compute: $Pr(S_2=x) = \frac{1}{2\pi}\int_{\mathbf{R}} e^{-itx}G_{S2}(\mathrm{e}^{it})dt$? $\endgroup$
    – Michel
    Commented Feb 10, 2013 at 14:12
  • $\begingroup$ Something in the way ... If we had a continuous distribution of which we can calculate the characteristic function (as you do), then this leads to a quick and nice result. In our case I need more time to think about it. There should be something easier. $\endgroup$
    – Richi W
    Commented Feb 10, 2013 at 16:13

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.