Unless you are using a symbolic language like Sage or Mathematica, you should always be aware that floating point numbers are… dangerous. In practice, the computer stores the digits and the position of the point. Which means that if a “decimal computer” (an imaginary computer that works with base 10 numbers) can store, say, 6 decimal digits, the number

` 1/3 = 0.333333`

is rounded down by about one millionth. But the operation:

` a = 1 + 1/3,000,000`

b = a - 1

leads to `b = 0.0000003`

, which means a rounding error of about 10%. Which makes you feel sorry you did not write:

` a = 1 - 1`

b = a + 1/3,000,000

Normally, this is not a big deal. But in probability computation, it is. Take the way the disjunction of probability is computed. Given `N`

events whose probability is `p`

, the probability that at least one of these events happens is:_{i}

**[1] **P(p_{1} or p_{2} or ... or p_{N}) = 1 - (1-p_{1})·(1-p_{2})·...·(1-p_{N})

It is immediately clear that if all `p`

are very small, and will be rounded like the 1/3 above, the error can become soon big enough to spoil your computation._{i}

And this is actually a pitty! Because:

**[2]** 1 - (1-p_{1})·(1-p_{2})·...·(1-p_{N}) = \Sum p_{i} + \Sum p_{i}·p_{j} - \Sum p_{i}·p_{j}·p_{k} ...

Where all sums run over 1 to N with indexes different among themselves. If we solve the formula algebraically, subtract 1 and only after substitute the values in the various p_{i}, the error would be smaller.

Let’s make an extreme example. If we had two events with p_{1}=p_{2}=10^{-6}:

**[3] **P(p_{1} or p_{2}) = 2·10^{-6} - 10^{-12} = 1.99999·10^{-6}

Imagine you have a “decimal computer” that can store only 3 digits plus the position of the point. If you solve **[3]** using the formula in **[2]**, even with such lousy computer the result would be a decent 1.99·10^{-6 }. This because both 2·10^{-6} and 10^{-12} do not need more than 3 digits. On the contrary, a naive computation, where the products are computed through **[1]** substituting the value of the variables, would give:

**[1'] **1 - 0.999·0.999 = 1.99·10^{-3}

A terribly wrong result. The short script below, in Python, solve the equation algebraically:

def probabilities_disjunction(probabilities: list): n = len(probabilities) # tree is a list of strings whose length is len(probabilities) with all possible 0/1 without the 0s. # e.g. for three probabilities: 001, 010, 011, 100, 101, 110, 111 tree = [] for i in range(1, 2**n): tree.append(str(bin(i))[2:].zfill(n)) v = [] for p in probabilities: v.append([1., -p]) prod = 0 # multiply all probabilities*(-1) according to tree for branch in tree: inter_prod = 1. for i, leaf in enumerate(branch): inter_prod *= v[i][int(leaf)] prod += inter_prod return -prod

It is clearly non optimised, and the number of computations grows exponentially with the number of probabilities.