[SciPy-dev] Sad sad sad... Was: Warning about remaining issues in stats.distributions ?

[Yaroslav Halchenko]

heh heh... very sad to see that the warning was simply ignored and 0.7.0
still has this issue on exactly the same command (hence advice to include
it to unittests was ignored as well):

>>> print scipy.__version__
0.7.0
>>> scipy.stats.rdist(1.32, 0, 1).cdf(-1.0+numpy.finfo(float).eps)
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/usr/lib/python2.5/site-packages/scipy/stats/distributions.py", line 117, in cdf
    return self.dist.cdf(x,*self.args,**self.kwds)
  File "/usr/lib/python2.5/site-packages/scipy/stats/distributions.py", line 625, in cdf
    place(output,cond,self._cdf(*goodargs))
  File "/usr/lib/python2.5/site-packages/scipy/stats/distributions.py", line 528, in _cdf
    return self.veccdf(x,*args)
  File "/usr/lib/python2.5/site-packages/numpy/lib/function_base.py", line 1886, in __call__
    _res = array(self.ufunc(*newargs),copy=False,
  File "/usr/lib/python2.5/site-packages/scipy/stats/distributions.py", line 525, in _cdf_single_call
    return scipy.integrate.quad(self._pdf, self.a, x, args=args)[0]
  File "/usr/lib/python2.5/site-packages/scipy/integrate/quadpack.py", line 185, in quad
    retval = _quad(func,a,b,args,full_output,epsabs,epsrel,limit,points)
  File "/usr/lib/python2.5/site-packages/scipy/integrate/quadpack.py", line 249, in _quad
    return _quadpack._qagse(func,a,b,args,full_output,epsabs,epsrel,limit)
  File "/usr/lib/python2.5/site-packages/scipy/stats/distributions.py", line 3046, in _pdf
    return pow((1.0-x*x),c/2.0-1) / special.beta(0.5,c/2.0)
ZeroDivisionError: 0.0 cannot be raised to a negative power

and my workaround doesn't work any more so I need to look for another
one.


On Tue, 09 Dec 2008, Yaroslav Halchenko wrote:

> > * distributions that have problems for some range of parameters
> so a good (imho) piece to add to unittests for the 'issues' to be fixed:

> scipy.stats.rdist(1.32, 0, 1).cdf(-1.0+numpy.finfo(float).eps)

> (tried on the SVN trunk to verify that it fails... discover the reason on
> your own ;-))

> For myself I resolved it with

>     __eps = N.sqrt(N.finfo(float).eps)
>     rdist = rdist_gen(a=-1.0+__eps, b=1.0-__eps, ....

> but I am not sure if that is the cleanest way... and may be some other
> distributions would need such tweakery to make them more stable.
-- 
                                  .-.
=------------------------------   /v\  ----------------------------=
Keep in touch                    // \\     (yoh@|www.)onerussian.com
Yaroslav Halchenko              /(   )\               ICQ#: 60653192
                   Linux User    ^^-^^    [175555]