Computing convolution density function with Maple

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

-
The PDF of the sum of random variables will be the convolution of the individual PDF's. Hence, take the Fourier Transform (FT) of the individual PDF's, multiply the FT's, and take the inverse FT of the product to obtain the convolution. – Rod Carvalho Sep 16 '12 at 18:26
I want to generalize that situation... I want to be able to change de distribution easily... I'm sure there is a way to do it like I want but I don't know how... – Jean-Philippe Le Cavalier Sep 16 '12 at 18:32

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]

-
Thanks for your answer. In facts, I have a convolution of 500 negative binomial distributions... I don't use Mathematica, there is probably a way to do it with Maple... – Jean-Philippe Le Cavalier Sep 16 '12 at 19:31