Computing convolution density function with Maple

Question

I would like to know how to compute a probability function of a convolution of Negative Binomial distribution with Maple.

Here is an easy example of what I want to do : '

With(Statistics):

X[1]:=RandomVariable(NegativeBinomial(2,0.5)):
X[2]:=RandomVariable(NegativeBinomial(6,0.3)):

S:=X[1]+X[2]:

ProbabilityFunction(S,0);

If I ask for the Mean, it works fine with Mean(S) but when I ask for the Probability Function, Maple gives me "FAIL" as answer.

Thank you!

Jean-Philippe

Answer 1

I am not really sure how to compute this with Maple, but here is how to get it in closed form.

Recall that the negative binomial distribution with parameters $n$ and $p$ has the following simple probability generating function: $$ \mathcal{P}_X\left(z\right) = \left( \frac{p}{1-(1-p)z} \right)^n $$ And that the probability generating function for the sum of two independent random variables $X_1+X_2$ is a product of individual generating functions: $$ \mathcal{P}_{X_1+X_2}(z) = \left( \frac{p_1}{1-(1-p_1)z} \right)^{n_1} \left( \frac{p_2}{1-(1-p_2)z} \right)^{n_2} $$ The probability mass function for $X_1+X_2$ is then read off as a series coefficient of PGF: $$ \mathbb{P}\left(X_1+X_2 = n\right) = [z^n] \mathcal{P}_{X_1+X_2}(z) $$ In particular: $$ \mathbb{P}\left(X_1+X_2=0\right) = p_1^{n_1} p_2^{n_2} = \mathbb{P}\left(X_1=0\right) \mathbb{P}\left(X_2=0\right) $$

By the way, here is what Mathematica gives for your problem:

In[14]:= Simplify[
 PDF[TransformedDistribution[
   x1 + x2, {Distributed[x1, NegativeBinomialDistribution[2, 1/2]], 
         Distributed[x2, NegativeBinomialDistribution[6, 3/10]]}], k]]

Out[14]= Piecewise[{{243*2^(-14 - k)*5^(-7 - k)*(66*5^(7 + k) - 
      15030*7^(3 + k) + 3*5^(7 + k)*k + 30673*7^(2 + k)*k - 
      4300*7^(2 + k)*k^2 + 60*7^(3 + k)*k^3 - 
            20*7^(2 + k)*k^4 + 2*7^(2 + k)*k^5), k >= 0}}, 0]

Compare with the series coefficient:

In[18]:= With[{p1 = 1/2, p2 = 3/10}, 
 Simplify[SeriesCoefficient[(p1^2*
      p2^6)/((1 + (-1 + p1)*z)^2*(1 + (-1 + p2)*z)^6), {z, 0, k}]]]

Out[18]= Piecewise[{{243*2^(-14 - k)*5^(-7 - k)*(66*5^(7 + k) - 
      15030*7^(3 + k) + 3*5^(7 + k)*k + 30673*7^(2 + k)*k - 
      4300*7^(2 + k)*k^2 + 60*7^(3 + k)*k^3 - 
            20*7^(2 + k)*k^4 + 2*7^(2 + k)*k^5), k >= 0}}, 0]