Numerical precision of product of probabilities (normal CDF)

Question

I'm trying to calculate $\prod_k{p_k}$ where $p_k$ are (potentially) very high probabilities of independent, zero-mean, standard normal random variables and $k>100$. However, I'm running into numerical problems using MATLAB (although the same problem occurs in Python/Scipy). Let's say $x=30$, then

normcdf(x)

returns 1, which is not precise. However, if I use normcdf(-x) (or normcdf(x,'upper')) instead, I get a value of 4.906713927148764e-198. I was hoping that I could then take 1 minus this value to get a more accurate probability. Unfortunately, the result gets rounded, as soon as I apply the subtraction:

>> normcdf(-x)
ans =
    4.906713927148764e-198
>> 1-(1-normcdf(-x))
ans =
     0

Is there any way to work around this issue?

Answer 1

If you really need to compute such tiny probabilities, one way is to use symbolic math and/or variable precision arithmetic.

For example, using the vpa function in the Symbolic Math toolbox:

X = sym(-300);
P = normcdf(X,0,1)
P2 = vpa(P)

which returns

P =

erfc(150*2^(1/2))/2


P2 =

7.449006262775352900552391145102e-19547

Of course converting this result back to floating-point (double(P2)) results in zero as P2 is less than eps(realmin). However, it's possible that if you do your calculations in variable precision and convert back to floating-point at the end you may be able gain a bit more accuracy. Just check to make sure that you're not wasting compute cycles.

Answer 2

In general, in 32-bit floating point arithmetic, $$1+\epsilon = 1,$$

where $|\epsilon| < |\epsilon_{\textrm{mach}}|$, the right-hand side being the machine epsilon.

I also question the need to compute a probability to 198 decimal places. There are about $10^{80}$ atoms in the universe. I cannot possibly imagine the need to compute probabilities so precisely that they could be used to predict something at the atomic level in two standard and one somewhat smallish universe.