1
$\begingroup$

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

  • 0
    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...2012-09-16

1 Answers 1

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] 
  • 0
    Thanks fo$r$ your answer. In facts, I have a convolution of 500 negative binomial $d$istributions... I $d$on't use Mathematica, there is probably a way to do it with Maple...2012-09-16