Function to compute the probability of exactly one event occurring out of N independent events

Question

Is there a built in function (or a function in one of the standard packages) that allows you to compute the probability of exactly one event occurring out of some known set of probabilities for N independent events?

Ideally what I'd want is something with a signature like:

Pr[independentProbabilities__] := ... 

To give an example, suppose I have a process that produces a letter from the set {A, B, C, D, X} with probabilities:

$$P(A) = 0.1\text{, } P(B) = 0.22\text{, } P(C) = 0.17\text{, } P(D) = 0.28$$

what is the probability that exactly one of the possibilities (A, B, C, or D) occurs?

Answer 1

Pr[independentProbabilities__] := 
 Block[{x, len = Length[{independentProbabilities}]},
  Probability[Sum[x[i], {i, len}] == 1, 
   Thread[x /@ Range[len] \[Distributed] 
     BernoulliDistribution /@ {independentProbabilities}]]
  ]

So

In[19]:= Pr[0.1, 0.22, 0.17, 0.28]

Out[19]= 0.414007

Or, doing the math

PrV2[independentProbabilities__] := 
 Total[Times @@@ ((1 - 
       ConstantArray[{independentProbabilities}, 
        Length[{independentProbabilities}]]) + 
     DiagonalMatrix[2 {independentProbabilities} - 1])]

Answer 2

This reply contributes a theoretical analysis of the question and provides a Mathematica solution that scales well for large lists of probabilities.

(Questions have appeared on our sister site, stats.stackexchange.com, involving arrays of up to a million probabilities, so this is of more than just academic interest.)

Theoretical solution

In clarifying comments to the question we learn that the argument independentProbabilities is a list of chances of independent events. By definition, when events $A$ and $B$ are independent, their probabilities multiply: $\Pr(A \cap B) = \Pr(A)\Pr(B)$. Moreover, (a) $A$ is also independent of the complement of $B$ ($B$ not happening) and (b) the probability of the complement of $B$ equals $1-\Pr(B)$. Finally, when events $E$ and $F$ cannot both happen, they are mutually exclusive and, to find the chance that at least one of them happens, their probabilities add: $\Pr(E \cup F) = \Pr(E) + \Pr(F)$.

From these (axiomatic) properties of probability we conclude that the chance of exactly one of independent events $A, B, \ldots, Z$ happening can be obtained by finding the chance that only $A$ happens, adding to that the chance that only $B$ happens, ..., etc.

The formula therefore is of the form

$$\eqalign{ \Pr(\text{Exactly one of }A,B,C,\ldots,Y,Z) &= \Pr(A)(1-\Pr(B))\cdots(1-\Pr(Z)) \\ &+ (1-\Pr(A))\Pr(B)(1-\Pr(C))\cdots(1-\Pr(Z)) \\ &+ (1-\Pr(A))\cdots(1-\Pr(Y))\Pr(Z). }$$

Mathematica implementation

An elegant way to compute this from an array of these probabilities exploits polynomial multiplication:

ClearAll[pr];
pr[p_List] := Block[{u, x = 1 - p, n = Length[p]},
  Coefficient[Times @@ (u + x), u] - n Times @@ x 
  ]

This works correctly even when an empty list is presented (the probability that exactly one of no events occurs should be 0).

Efficiency is an issue with some solutions

E.g.,

data = RandomReal[{0, 1}, {100000, 4}];
Timing[pr /@ data;]

takes 5.0 seconds. Compare this to 94 seconds for a solution involving BernoulliDistribution. Changing the test data from 100,000 lists of 4 probabilities to 1000 lists of 400 probabilities increases the time to 12.4 seconds, but the solution based on BernoulliDistribution is unable to finish! Extrapolating its performance for smaller problems suggests it would need about 2000 years.

Analysis of alternatives

The matrix-based PrV2 solution offered by @rojo is quite efficient--it works about twice as quickly as the solution offered here does--until you try passing longish arrays of probabilities as arguments. The break-even point is around 120 probabilities. After this, PrV2 scales quadratically in the argument size $n$, whereas PrV2 scales as $O(n\log(n))$. This will be a problem with any direct implementation of the formula, because it involves a sum of $n$ products of $n$ terms each: obviously $O(n^2)$ effort. The polynomial-based algorithm is superior for large lists of probabilities because (apparently) the polynomial algorithms use efficient convolutions, which tend to be $O(n\log(n))$ operations rather than $O(n^2)$.

For extremely fast calculation, you can obtain $O(n)$ asymptotics provided none of the arguments equals $0$ (or is so small as to cause numerical overflow):

Clear[pr3];
pr3[p_List] := Block[{x}, (Times @@ (1 - p)) Sum[x/(1 - x), {x, p}]]

Under these slightly limited conditions, this solution will be superior in execution time to any of the others, no matter what the length of its argument:

data = RandomReal[{0, 1}, {100000, 4}];
Timing[pr3 /@ data;]

(1.7 seconds)

data = RandomReal[{0, 1}, {1, 400000}];
Timing[pr3 /@ data;]

(0.13 seconds)