The weighted sum of two independent Poisson random variables

Question

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?

Answer 1

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.

Answer 2

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$

Answer 3

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.