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.
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