Avoiding underflow for joint probabilities using NumPy

Question

I face a problem of estimating a joint probability for independent variables in a simple setup. Currently I have an array of 100 random variables and I would like to obtain their joint probability without failing into the underflow problem. Any ideas how to achieve this goal in numpy? If possible?

If not could someone please explain me further the role of NumPy routine (logaddexp) as I thought it might be a help for me in such situation.

Answer 1

logaddexp lets you expand the range (reducing the precision) of representable values by storing and processing their logarithm instead.

e1, e2 = np.log(p1), np.log(p2)    # convert p1, p2 to log form
e1e2 = e1 + e2                     # e1e2 is now np.log(p1 * p2)
e3 = np.logaddexp(e1, e2)          # e3 = np.log(p1 + p2)

You just need to go through your code changing ** to *, * to + and + to np.logaddexp, and convert back with np.exp at the end.

The normal 64-bit double-precision floating point has least positive normal value 2.2E-308; storing logs gives you an effective least positive normal 1E-(1.7E308).