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

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?

-
I don't really understand the question. If you draw one letter, then of course exactly one will occur. Also, those probabilities you gave don't sum up to 1, is the remaining associated with X? – Szabolcs Mar 15 '12 at 7:16
One possible interpretation is that $N$ independent draws are made from the given distribution. The probability that exactly one of $A, \ldots, D$ occurs is the chance that $X$ occurs exactly $N-1$ times: that's given by the binomial probability With[{p = Plus @@ independentProbabilities}, N p(1-p)^(N-1)]. But this is so simple I wonder whether it is the intended interpretation... – whuber Mar 15 '12 at 16:16
@whuber - In this particular case, you are correct - but that is more a flaw in my example. What I'm looking for is a functional expansion of the identity: P(A ^ B ^ C) == (P(A)*(1-P(B))*(1-P(C)) + ((1-P(A))*P(B)*(1-P(C))) + ((1-P(A))*(1-P(B))*P(C)). – LBushkin Mar 16 '12 at 5:50
If by "^" you mean conjunction, then that identity appears to be incorrect or meaningless. I doubt it refers to a single outcome of the process, because then "A^B^C" is the empty set (since outcomes A, B, and C are mutually exclusive). If it refers to a sequence of three independent outcomes, we could interpret the lhs as the probability of the sequence ABC, but the rhs gives the sum of the chances that A occurs first or B second or C third; it is going to be greater than the lhs, because it includes outcomes of the form AXX, XBX, XXC, ABX, etc. So what did you really mean to write? – whuber Mar 16 '12 at 13:59
@whuber: Let me see if I can put it a different way. If I told that you that there is a 10% chance that it will rain tomorrow, a 22% chance that there will be a sale at Macy's, and a 17% chance that your flight will be delayed, could you tell me the probability of exactly one and only one of those things occurring? – LBushkin Mar 17 '12 at 5:59

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])]

-

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)

-
Great stuff. You'd do well at Project Euler, if you're not already there. – Mr.Wizard Mar 19 '12 at 22:23
Incidentally you don't need Block in the definition of pr3. – Mr.Wizard Mar 19 '12 at 22:25
Mr. Wizard, I believe you, but frequently I run into trouble when I don't protect index variables in Sum, Table, or whatever, even though the docs promise that there is an implicit block already there. That's probably due to some misunderstanding of mine, but in the meantime discretion seems to be the better part of valor... – whuber Mar 19 '12 at 22:41
I would like to see an example of these running into trouble. – Mr.Wizard Mar 19 '12 at 23:02
I knew you would respond with that :-). I'll keep this site in mind the next time I run into apparent difficulties with Mathematica's scoping and I'll start a new question rather than continuing this comment thread. – whuber Mar 20 '12 at 15:40