dist_math.incomplete_beta is very slow

[ricardoV94]

I notice this when running test_distributions.py. Any distribution logcdf method that uses incomplete_beta takes ages to run.

Right now it is used in the logcdf methods of Beta, StudentT, Binomial, NegativeBinomial, and ZeroInflated versions of the last two distributions.

Is there a reason why this cannot be implemented as a Theano C Op? Scipy has C implementation here, which seems to be doing exactly the same as our algorithm in dist_math.py:
https://github.com/scipy/scipy/blob/master/scipy/special/cephes/incbet.c

In addition, this would probably make it trivial to implement a vectorized version for tensors, allowing the logcdf methods of the distributions above to evaluate more than one value at a time.

[ricardoV94]

Looking at this #4454 CI run: https://github.com/pymc-devs/pymc3/pull/4454/checks?check_run_id=1807188616

The top 5 slowest distribution tests include all of the ones mentioned above (with the notable exception of the StudentT) and together took nearly 1 hour to complete.

[ricardoV94]

The StutentT is not in that list because the number of values tested by the logcdf was already reduced to 10 before:

https://github.com/pymc-devs/pymc3/blob/89916ad724ae88ebb4c1af3d65919b87b318f180/pymc3/tests/test_distributions.py#L1131-L1137

[junpenglao]

Is the gradient available? Looks like it might be much faster create a theano Op to call scipy directly similar to
https://github.com/pymc-devs/pymc3/blob/3cd9eb3004956f3624f187ff428b8e33dc4aa729/pymc3/distributions/dist_math.py#L369-L382

[ricardoV94]

I am not sure about the gradient, I will do some research. If anyone knows please chime in :)

[ricardoV94]

I found some formulae for the partial derivatives here: https://functions.wolfram.com/GammaBetaErf/Beta3/20/01/

The partial derivative with respect to x is trivial (given the beta is defined as an integral up to x)

The partial derivatives with respect to a, and b are more tricky. They require the Hypergeometric function 3F2, which is not given by scipy AFAIK. It seems to exist in mpmath though: https://mpmath.org/doc/current/functions/hypergeometric.html#hyp3f2.