Does Mathematica have an implementation of the Poisson binomial distribution?

Question

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?

Answer 1

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 ]

Answer 2

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...

Answer 3

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}$

Answer 4

I have answered this question using a full implementation of the Poisson Binomial.