[SciPy-User] Question on scipy.special.gammaincc

[Skipper Seabold]

On Wed, May 25, 2011 at 5:28 AM, Pauli Virtanen <pav at iki.fi> wrote:
> Tue, 24 May 2011 21:59:29 -0400, Skipper Seabold wrote:
> [clip]
>> Indeed. For the time being, would a patch to return nan be useful? I
>> don't see how to correctly raise exceptions from C, unless someone has a
>> pointer to an example.
>
> A patch would be useful. In addition to returning `nan`, call
>
>        mtherr("gammaincc", DOMAIN)
>
> from "cephes/mconf.h", which will raise a warning (and maybe an error)
> depending on what the user has set up.
>

Right, thanks. I was also wondering if it's possible to raise
something like a NotImplementedError for negative real numbers since
only the negative integers and zero are poles. I did as you suggested
though.

I also think the cephes docs are slightly misleading for incomplete
gamma. It says it only works for positive arguments for both, but I
think the integration limit should be allowed to be zero (the
chi-square and gamma distribution functions depend on this) no matter
what the other arguments are. I also updated the docs.

Can someone check this small patch?

https://github.com/jseabold/scipy/compare/master...gammainc-nanfix

Versus mpmath:

[20]: lower_int = 0

[cephes]
[21]: special.gammainc(0,1)
[21]: nan

[cephes]
[22]: mpmath.gammainc(0,lower_int,1, regularized=True)
[22]: mpf('+inf')

[cephes]
[23]: special.gammainc(1,0)
[23]: 0.0

[cephes]
[24]: mpmath.gammainc(1,lower_int,0, regularized=True)
[24]: mpf('0.0')

[cephes]
[25]: special.gammainc(0,0)
[25]: 0.0

[cephes]
[26]: mpmath.gammainc(0,lower_int,0, regularized=True)
[26]: mpf('0.0')

[cephes]
[27]: special.gammainc(1,1)
[27]: 0.63212055882855778

[cephes]
[28]: mpmath.gammainc(1,lower_int,1, regularized=True)
[28]: mpf('0.63212055882855767')


Skipper