numpyscipyprecisionfloating-accuracympmath

How to handle very small numbers (<1e-15) in Python, passed to Scipy/Numpy functions?


I would like to compute some scipy functions, such as: scipy.special.gammaincinv (https://docs.scipy.org/doc/scipy/reference/generated/scipy.special.gammaincinv.html#scipy.special.gammaincinv).

The challenge, however, is that one of the arguments I need to use is 1 - a where a is a very small number, a=1e-18 in this case.

Therefore, as Python's small number precision is limited to about 1e-15, Python computes 1 - a = 1 - 1e-18 = 1.0, i.e. it's inaccurate.

I could use libraries like mpmath to get an accurate calculation of 1 - a but then the returned object is not a normal number, so can't be passed to numpy / scipy functions.

Therefore, my question is: how best to use scipy & numpy functions when the input numbers include precision beyond that of standard floats (ideally, not simply rewriting the functions themselves using mpmath-type constructs?)

EDIT: The exact expression I want to calculate is: gammaincinv(10, 1-a) where a is a very small number (a=1e-18). Due to limited precision, 1-a is calculated as 1.0, which gives the wrong value.


Solution

  • Use gammainccinv(10, a) instead of gammaincinv(10, 1 - a). This is exactly what gammainccinv is for.

    For moderate values of a, you can see that they give the same value:

    In [14]: from scipy.special import gammaincinv, gammainccinv
    
    In [15]: a = 0.125
    
    In [16]: gammaincinv(10, 1 - a)
    Out[16]: 13.688247592108494
    
    In [17]: gammainccinv(10, a)
    Out[17]: 13.688247592108494
    
    
    

    For small a, all precision is lost in the subtraction 1 - a

    In [18]: a = 1e-18
    
    In [19]: gammaincinv(10, 1 - a)
    Out[19]: inf
    

    Using gammainccinv avoids that subtraction and gives the correct answer:

    In [20]: gammainccinv(10, a)
    Out[20]: 66.57188367091577