# Does Mathematica have an implementation of the Poisson binomial distribution?

I need to work out the probability of having $$k$$ successful trials out of a total of $$n$$ when success probabilities are heterogeneous. This calculation relates to the Poisson Binomial Distribution. Does Mathematica, or perhaps the Mathstatica add-on, have an implementation for that?

Mathematica does not know about the PoissonBinomialDistribution, but you can use the formula given for the PDF on Wikipedia:

PoissonBinomialDistribution[ plist : { __?NumericQ } ] := With[
{
n = Length @ plist,
c = Exp[(2 I \[Pi])/(Length@plist + 1)]
}
,
ProbabilityDistribution[
Re[ 1/(n + 1) Sum[c^(-l k) Product[1 + (c^l - 1) plist[[m]] , {m, 1, n }], {l, 0, n}] ]
,
{k, 0, n, 1}
]
] /; AllTrue[ plist, 0 <= # <= 1& ]


Now we may model a quality control where fault type 1 has a prob of 4% and fault types 2 and 3 have a prob of 7%:

dist = PoissonBinomialDistribution[ {0.04, 0.07, 0.07} ];


With this we find the probability for 3 faults:

Probability[ k == 3, k \[Distributed] dist ]// PercentForm


0.0196 %

Update

Another possibility is to work with TransformedDistribution. Doing so also allows for symbolic evaluation:

PoissonBinomialDistribution::fargs = "Probabilities must be between zero and one";
PoissonBinomialDistribution[ plist : {__?NumericQ | __Symbol } ] := With[
{
xList = Table[ Unique["x"], Length @ plist ]
}
,
If[ Not @ AllTrue[ plist, 0 <= # <= 1 &],
Return[ Message[PoissonBinomialDistribution::fargs]; \$Failed ]
];
TransformedDistribution[
Total @ xList,
Thread[ xList \[Distributed] Map[ BernoulliDistribution, plist] ]
]
]

PDF[ PoissonBinomialDistribution[{p1, p2, p3}], x ]


• This produces wildly incorrect results pretty quickly... – ciao Apr 25 at 5:59
• @ciao Will the TransformedDistribution solution I just added also produce incorrect results? – gwr Apr 25 at 6:02
• The latter s/b fine, but will get very slow with more than a couple dozen probabilities. – ciao Apr 25 at 6:38
• @ciao Indeed, I just tested with a bunch of random numbers: Already with about 15 probs it takes "hours". :-( – gwr Apr 25 at 6:41
• @gwr your original solution is much faster. – user120911 Apr 25 at 7:11

From what it appears you require, the following should suffice:

pb = Fold[ListConvolve[##, {1, -1}, 0] &, Transpose[{1 - #, #}]] &;


It should provide more accurate results for reasonably sized probability vectors.

As an example, a vector of 100 probabilities:

SeedRandom@1234
testps = RandomReal[{0, 1}, 100];


Timings:

r1 = pb@testps; // AbsoluteTiming // First

r2 = Re@PDF[PoissonBinomialDistribution[testps], Range[0, Length@testps]]; // AbsoluteTiming // First


0.0008377

0.137304

Comparing accuracy with a random sample of the PMFs (pb left, PoissonBinomialDistribution right. The correctness of pb was verified for full PMF):

With[{p = RandomSample[Range@100, 10]},
Transpose[{r1[[p]], r2[[p]]}]] // TableForm


You can of course use the same to generate generic results for n probabilities, à la Carl Woll's answer, and then replace probabilities as desired:

np = 18;

r1 = Table[p[np, x], {x, 0, np}]; // AbsoluteTiming // First

r2 = pb[Array[p, np]]; // AbsoluteTiming // First

%%/% // Round

And @@ MapThread[FullSimplify[#1 == #2] &, {r1, r2}]


504.622

0.0738808

6830

True

Same results, but ~1/7000th the time taken. I tried with np of 20, gave up waiting, probably 1/20000th the time or so...

• (+1) Very nice. If one uses machine precision numbers and adds Chop and Re for the Fourier-Transform-Solution the timing can at least be brought down to around 0.034 secs on my machine (which is of course still way off this solution). – gwr Apr 25 at 10:22
• Even faster should be Threshold instead of Chop... – Henrik Schumacher Apr 29 at 20:15

Another possibility is:

p[n_, k_] := SeriesCoefficient[Product[1 - (1 - ϵ) p[i], {i, n}], {ϵ, 0, k}]


Reproducing gwr's results:

Table[p[3, i], {i, 0, 3}] //Simplify //Column //TeXForm


$$\begin{array}{l} -(p(1)-1) (p(2)-1) (p(3)-1) \\ -2 p(3) p(2)+p(2)+p(3)+p(1) (-2 p(3)+p(2) (3 p(3)-2)+1) \\ p(2) p(3)+p(1) (-3 p(3) p(2)+p(2)+p(3)) \\ p(1) p(2) p(3) \\ \end{array}$$