[Numpy-discussion] BigInteger equivalent in numpy

[josef.pktd at gmail.com]

On Sun, Jun 7, 2009 at 6:05 AM, David Cournapeau
<david at ar.media.kyoto-u.ac.jp> wrote:
> (Please do not send twice your message to numpy/scipy ML, thank you)
>
> wierob wrote:
>> Hi,
>>
>> int64 and float seem to work for the stderr calculation. Now, the
>> calculation of the p-value causes an underflow.
>>
>> File "C:\Python26\lib\site-packages\scipy\stats\distributions.py", line 2829,in _cdf
>>     return special.stdtr(df, x)
>> FloatingPointError: underflow encountered in stdtr
>>
>> It seems that the error occurs in scipy.special.stdtr(df, x) if df =
>> array([13412]) and x = array([61.88071696]).
>>
>
> The error seems to happen in the cephes library (the library we use for
> many functions in scipy.special, including this one). I don't know
> exactly what triggers it (your df is quite high, though, and the cdf is
> already quite past the 1-eps at your point).
>

The result is correct, and we don't get the exception at the default
warning/error level, so this looks ok to me.

>>> stats.t._cdf(61.88071696, 13412)
1.0
>>> stats.norm._cdf(61.88071696)
1.0
>>> stats.norm._sf(61.88071696)
0.0
>>> stats.norm._ppf(1-1e-16)
8.2095361516013874

Using a double precision library won't help it find out whether the
answer is (1 - 1-e30) or (1 - 1-e80)

I made comments on a similar case yesterday where the overflow
actually shows up in the result without any warning or exception.


http://projects.scipy.org/scipy/ticket/962

But I only know about floating point precision what I learned in these
news groups, and playing with limiting cases in stats.distributions.

using error level to raise doesn't seem useful if you want for example
to work with include floating point inf

>>> np.seterr(all="raise")
{'over': 'ignore', 'divide': 'ignore', 'invalid': 'ignore', 'under': 'ignore'}
>>> x = np.array([0,1])
>>> x
array([0, 1])
>>> x/x[0]
Traceback (most recent call last):
  File "<pyshell#140>", line 1, in <module>
    x/x[0]
FloatingPointError: divide by zero encountered in divide


>>> x = np.array([0.,1.])
>>> x/x[0]
Traceback (most recent call last):
  File "<pyshell#143>", line 1, in <module>
    x/x[0]
FloatingPointError: divide by zero encountered in divide
>>> np.seterr(all="ignore")
{'over': 'raise', 'divide': 'raise', 'invalid': 'raise', 'under': 'raise'}
>>> x/x[0]
array([ NaN,  Inf])

Josef