Timeline for How can I efficiently model the sum of Bernoulli random variables?
Current License: CC BY-SA 3.0
4 events
| when toggle format | what | by | license | comment | |
|---|---|---|---|---|---|
| Mar 24, 2012 at 18:08 | comment | added | whuber♦ |
Here's the Mathematica code: multinomial[p_] := Module[{lc, condense}, lc = Function[{s}, ListConvolve[s[[1]], s[[2]], {1, -1}, 0]]; condense = Function[{s}, Map[lc, Partition[s, 2, 2, {1, 1}, {{1}}]]]; Flatten[NestWhile[condense, Transpose[{1 - p, p}], Length[#] > 1 &]]] To apply it, do something like p = RandomReal[{0, 1}, 40000]; pp = multinomial[p];. This creates the probabilities p and then computes the exact distribution pp. NB When the mean of p is not extreme, the distribution is very close to normal: that leads to a much faster algorithm yet.
|
|
| Mar 24, 2012 at 18:05 | comment | added | whuber♦ | (+1) I found a way to make this idea work. It requires two techniques: first, use FFT for the convolutions; second, don't do them sequentially, but divide and conquer: do them in disjoint pairs, then do the results in disjoint pairs, etc. The algorithm now scales as $O(n\log n)$ rather than $O(n^2)$ for $n$ probabilities. For instance, Mathematica can compute the entire distribution for 40,000 probabilities in just 0.4 seconds. (1,000,000 are calculated in 10.5 seconds.) I will provide code in a followup comment. | |
| Mar 23, 2012 at 19:17 | comment | added | whuber♦ | Have you tried this with 40K variables?? (I wonder how many hours or days of computation it takes...) | |
| Mar 23, 2012 at 17:50 | history | answered | alex | CC BY-SA 3.0 |