Multiple discrete RV convolutions performance

Question

I have a discrete random variable with PDF of (just a simplified example)

Piecewise[{{1/3, x == 5 || x == 10 || x == 15}}, 0]

Actually, in my code this is a MixtureDistribution of a number of discrete distributions.

What I need is to obtain a multiple convolution (for example a 100-fold) of independent identically distributed random variables, distributed according to the provided distribution.

I tried all known to me ways:

1) TransformedDistribution

TransformedDistribution[Plus @@ Flatten[Table[tX[v], {v, 0, vMax}]], Flatten[Table[tX[v] \[Distributed] resF, {v, 0, vMax}]]]

2) multiple DiscreteConvolve

DiscreteConvolve[PDF[resF, x], PDF[resF, x], x, y]

where resF is the provided initial distribution.

The first method provides the needed result, but it takes way too long and the time grows exponentially - approx 60 seconds to obtain a ~10-fold convolution, a few minutes for the next.

The second method seem to provide the correct answer for a couple of first convolutions (and it take longer compared to the first method), but the next convolutions give 0.

Just to compare I wrote a simple head-on program on C# to do the task and it takes only ~2-3 second to do the 300-fold convolution. And there is way to optimize it by removing values, higher than the threshold on the go.

Ultimately, I need the probability that the final convolution will be smaller than a provided threshold, i.e. CDF.

Can someone tell me why is it taking so long or if I'm doing this the wrong way.

Thanks!

EDIT

Some code that I test:

ClearAll[resFF, vMax];
AbsoluteTiming[
 resF = EmpiricalDistribution[{1/3, 1/3, 1/3} -> {5, 10, 15}];
 vMax = 10;
 PDF[TransformedDistribution[
   Plus @@ Flatten[Table[tX[v], {v, 1, vMax}]], 
   Flatten[Table[tX[v] \[Distributed] resF, {v, 1, vMax}]]], x]
]

Here resF is the distribution I mentioned in the beginning, and vMax is the number of convolutions I want to obtain. So this code gives the PDF of a 10-fold convolution.

On my machine it takes 3 seconds to do the 10-fold convolution, 8 seconds for 11-fold and 25 seconds for 12-fold. As I need more than 100-fold one, this performance is not sufficient.

EDIT 2

As Jim Baldwin said, I missed the fact that I can obtain the result using the generating functions of discrete random variables.

Thus, the code I came up with that gives me the desired value of vMax-fold convolution CDF in point 1000 is:

ClearAll[resF, gen, vMax];
vMax = 150;
resF[x_] := Piecewise[{{1/3, x == 5 || x == 10 || x == 15}}, 0];
gen = GeneratingFunction[resF[n], n, x];

Total@Select[CoefficientList[(gen)^vMax, x][[1 ;; 1000]], # > 0 &]

However, in this case I have to scale the initial data, if the outcomes of RV are not integer.

PS: This is out of my own curiosity, but still - if anyone know why the TransformedDistribution solution works so slow, please, tell me. Thanks!

Answer 1

Probability generating functions might be the way to go. Here's the pgf (probability generating function) of a single random variable:

pgf = t^5/3 + t^10/3 + t^15/3

The probability generating function for the sum of n independent and identically distributed random variables is

pgf^n

You can get the probability for any value of the resulting sum (as I assume what you mean by the convolution is that you want the distribution of the sum of the n random variables) by using CoefficientList and a bit of clean-up:

pgf = t^5/3 + t^10/3 + t^15/3;
cl = CoefficientList[pgf^10, t];
x = Range[0, Length[cl] - 1];
pdf = Transpose[{x, cl}];
Select[pdf, #[[2]] > 0 &]]

{{50, 1/59049}, {55, 10/59049}, {60, 55/59049}, {65, 70/19683}, {70, 
  205/19683}, {75, 484/19683}, {80, 950/19683}, {85, 1580/19683}, {90,
   2255/19683}, {95, 8350/59049}, {100, 8953/59049}, {105, 8350/
  59049}, {110, 2255/19683}, {115, 1580/19683}, {120, 950/
  19683}, {125, 484/19683}, {130, 205/19683}, {135, 70/19683}, {140, 
  55/59049}, {145, 10/59049}, {150, 1/59049}}

As far as timing...Here is an example for n=100:

Timing[pgf = t^5/3 + t^10/3 + t^15/3;
 cl = CoefficientList[pgf^100, t];
 x = Range[0, Length[cl] - 1];
 pdf = Transpose[{x, cl}];
 Select[pdf, #[[2]] > 0 &];]
(* {0.0312002`,Null} *)