16
$\begingroup$

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?

$\endgroup$
16
$\begingroup$

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 ]

Symbolic PDF

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

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

enter image description here

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

$\endgroup$
  • $\begingroup$ (+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). $\endgroup$ – gwr Apr 25 at 10:22
  • $\begingroup$ Even faster should be Threshold instead of Chop... $\endgroup$ – Henrik Schumacher Apr 29 at 20:15
3
$\begingroup$

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

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.