Mrthdrb

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 ] := 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 ]

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

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

$beginarrayl
-(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) \
endarray$