1
$\begingroup$

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?

$\endgroup$
11
  • $\begingroup$ What are the mean and variances? $\endgroup$
    – Samrat Mukhopadhyay
    May 23, 2014 at 15:53
  • $\begingroup$ @SamratMukhopadhyay The variables are normalized to zero mean and unit-variance (I updated the question) $\endgroup$
    – ws6079
    May 23, 2014 at 15:55
  • $\begingroup$ In what application is $1-4.9\times 10^{-198}$ a necessary level of precision??? $\endgroup$
    – Emily
    May 23, 2014 at 15:56
  • $\begingroup$ @Arkamis the idea is that the product of the $k$ probabilities will give me $P( \cap_k \{X \le x_k\})$ which should be in the order of $1 \times 10^{-6}$ $\endgroup$
    – ws6079
    May 23, 2014 at 16:02
  • 1
    $\begingroup$ @ws6079: I think you mean that they return zero. Try log1p(-normcdf(-x)), which, in this case, returns exactly -normcdf(-30). eps(1) is much larger than normcdf(-30) so you'll never be able to subtract it in floating point. $\endgroup$
    – horchler
    May 23, 2014 at 21:21

2 Answers 2

Reset to default
1
$\begingroup$

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.

$\endgroup$
1
  • $\begingroup$ Thanks - this is what I was looking for! $\endgroup$
    – ws6079
    May 23, 2014 at 22:21
0
$\begingroup$

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.

$\endgroup$
1
  • $\begingroup$ I'm trying to calculate the joint probability of independent events, each one of which has a high probability of occurring - well specifically I'm interested in the complement of this probability. I would be grateful for any suggestions for alternative approaches. $\endgroup$
    – ws6079
    May 23, 2014 at 16:15

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .