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!
rv
??) $\endgroup$ – Szabolcs Oct 28 '16 at 14:28tX[v]
, norvMax
, norresF
. It's just gobbledegook at the moment, leaving the reader to guess. $\endgroup$ – wolfies Oct 28 '16 at 15:00Coefficient
orCoefficientList
to get the probabilities of interest. $\endgroup$ – JimB Oct 28 '16 at 15:04