0

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.

Geoff2009max
  • 3
  • 2
  • The best way, when feasible, is to decompose your expression so that the 1 and `a` are processed separately. Thus your question needs to show exactly the expression you wish to calculate, not just "such as" – Reinderien Feb 21 '23 at 15:50
  • Thanks. This is a good point. I tried to edit the Q but SO has too many pending edits. 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. I wonder how I might implement your suggestion of processing these separately please? – Geoff2009max Feb 21 '23 at 16:34
  • `ufunc` like this use compiled code, so inputs are limited to the C types used in compilation (usually float and double). – hpaulj Feb 21 '23 at 16:40
  • Thanks. So this sounds like it can't be done? I'm sure others have come across similar problems with precision when wanting to use numpy/scipy functions so I wonder if there are any hacks or workarounds? – Geoff2009max Feb 21 '23 at 18:07
  • 2
    for values that I tested `sc.gammainccinv(10, x)` is the same as `sc.gammaincinv(10,1-x)` – hpaulj Feb 21 '23 at 21:01
  • Really? I get very different results, e.g.: `from scipy.special import gammaincinv` then `a = 1e-18` initially. Then `gammaincinv(10, 1-a)` returns `inf` while `gammaincinv(10, a)` returns `0.07224835788588542`. You can see the problem of finite precision by doing `gammaincinv(10, 0.99999)` then ``gammaincinv(10, 0.99999999999)` - note the output changes. But after 15 decimal points, Python rounds the 2nd argument to 1.0, which gives an erroneous `inf` result. – Geoff2009max Feb 22 '23 at 10:33

1 Answers1

0

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
Warren Weckesser
  • 110,654
  • 19
  • 194
  • 214