Return NaN in C implementations of SciPy Ops

[ricardoV94]

This is related to pymc-devs/pymc#4340

I am trying to use the GammaIncc Op, which makes use of the C code function GammaQ, which begins with an assertion for the two parameters. This is problematic when using the tensor version of the scalar op (via @_scal_elemwise), as I don't see how to gracefully avoid the assertion error.

To make this concrete, gammaincc is being used in the logcdf method of the pymc3 InverseGamma distribution. It is wrapped inside a bound which seems to lazy evaluate when working with scalars, but not with tensors. This allows it to gracefully return -np.inf when evaluated with invalid scalar arguments, but not when using a tensor/array.

The scipy.special.gammaincc(x, k) counterpart does not have this issue.

Does anyone know a workaround to avoid this issue when working with tensors?

pm.InverseGamma.dist(alpha=-1, beta=5).logcdf(1).eval()  # returns -np.inf
pm.InverseGamma.dist(alpha=-1, beta=5).logcdf(np.array([1])).eval()  # fails assertion and exits process

python3: /home/ricardo/.theano/compiledir_Linux-5.4--generic-x86_64-with-glibc2.29-x86_64-3.8.5-64/tmpf82kd0bh/mod.cpp:240: double GammaQ(double, double): Assertion `(n > 0) && (x >= 0)' failed.
Aborted (core dumped)

https://github.com/pymc-devs/pymc3/blob/043129243436b875a9ab18180549491c4678b7c6/pymc3/distributions/continuous.py#L2654-L2659

return bound(
    tt.log(tt.gammaincc(alpha, beta / value)),
    value >= 0,
    alpha > 0,
    beta > 0,
)

https://github.com/pymc-devs/Theano-PyMC/blob/2601e7acf8442422cab19db5bd029a5581f0bdec/theano/tensor/basic.py#L2529-L2531

https://github.com/pymc-devs/Theano-PyMC/blob/819122e66bcd089b2a3f472542800feead0bf755/theano/scalar/basic_scipy.py#L598-L639

https://github.com/pymc-devs/Theano-PyMC/blob/819122e66bcd089b2a3f472542800feead0bf755/theano/scalar/c_code/gamma.c#L227-L233

[dfm]

This is a good point! While it is true that the results are undefined for those parameters, I think that all the functions in gamma.c should be updated to return NaN for invalid parameters rather than a C-assertion.