Replace slow dist_math::incomplete_beta with aesara op

[ricardoV94]

This PR replaces the incomplete_beta method in dist_math by an aesara op betainc that wraps the scipy.special.betainc method.

This provides a large increase in speed (and graph compilation) even without providing the c_code implementation directly (this can always be added later). It also allows the evaluation of the logcdf for multiple values in Beta, StudentT, Binomial, NegativeBinomial (and zero inflated variants). I removed the previous limitation warning when trying to evaluate these methods for multiple values. These changes are critical for any practical implementation of Truncated Versions of these distributions.

Gradients were adapted from the STAN implementation, and require the use of two extra "pseudo aesara ops" that simply wrap the equivalent python functions. This does not allow for further graph optimizations, but I doubt aesara could do much with the previous scan implementation either. Let me know if you think there is a better solution.

Depending on what your PR does, here are a few things you might want to address in the description:

• what are the (breaking) changes that this PR makes? dist_math::incomplete_beta was deprecated in favor of the new betainc op
• important background, or details about the implementation
• are the changes—especially new features—covered by tests and docstrings?
• linting/style checks have been run
• consider adding/updating relevant example notebooks
• right before it's ready to merge, mention the PR in the RELEASE-NOTES.md

Closes #4420

[OriolAbril]

not sure if this should be merged into master or into v4 🤔 cc @brandonwillard @twiecki @michaelosthege

[brandonwillard]

It shouldn't be difficult to rebase v4, so we could merge something like this. I think the problem would have more to do with our release plan (e.g. are we going to make version 3 releases with the theano -> aesara change?)

[ricardoV94]

Strange failure on 32-bit run, for the logcdf of the beta.

E               AssertionError: 
E               Arrays are not almost equal to 3 decimals
E               {'alpha': array(100., dtype=float32), 'beta': array(0.1, dtype=float32), 'value': array(0.001, dtype=float32)}
E               x and y -inf location mismatch:
E                x: array(-inf, dtype=float32)
E                y: array(-697.172)

I don't understand why the aesara op would underflow but not the scipy one. It seems scipy is basically calling betainc and then taking the logarithm, which is exactly the same thing we are doing:

https://github.com/pymc-devs/pymc3/blob/53b642e4124feb0092ef5492685f44420f80cf3b/pymc3/distributions/continuous.py#L1319

https://github.com/scipy/scipy/blob/master/scipy/stats/_continuous_distns.py#L610-L611

Afaik the scipy should also underflow to 0 and then -inf in a float-32 environment:

np.log(st.beta(100.0, 0.1).cdf(0.001).astype('float32'))
<ipython-input-32-b9435d6250d7>:1: RuntimeWarning: divide by zero encountered in log
  np.log(st.beta(100.0, 0.1).cdf(0.001).astype('float32'))
Out[32]: -inf

[michaelosthege]

@ricardoV94 this error could indicate that y ran with float64:

E                x: array(-inf, dtype=float32)
E                y: array(-697.172)

[ricardoV94]

@ricardoV94 this error could indicate that y ran with float64:

E                x: array(-inf, dtype=float32)
E                y: array(-697.172)

I was probably unclear. The error appeared in the 32bit GH run: https://github.com/pymc-devs/pymc3/runs/2075254035?check_suite_focus=true

AFAIK both should have underflowed.

[ricardoV94]

In terms of V3 / V4, I think it's fine if this only appears in the V4 branch. I can toggle this as a draft until we are ready for it.

[ricardoV94]

Waiting a bit on this one following some discussion on the STAN side to make sure these gradient approximations are good enough: stan-dev/math#2417