Replace incomplete_beta with betainc and speedup logcdf tests

RELEASE-NOTES.md

@@ -3,12 +3,13 @@
 ## PyMC3 vNext (4.0.0)
 ### Breaking Changes
 - ⚠ Theano-PyMC has been replaced with Aesara, so all external references to `theano`, `tt`, and `pymc3.theanof` need to be replaced with `aesara`, `at`, and `pymc3.aesaraf` (see [4471](https://github.com/pymc-devs/pymc3/pull/4471)).
+- ⚠ PyMC3 now requires Scipy version `>= 1.4.1` (see [4736](https://github.com/pymc-devs/pymc3/pull/4736)).
 - ArviZ `plots` and `stats` *wrappers* were removed. The functions are now just available by their original names (see [#4549](https://github.com/pymc-devs/pymc3/pull/4471) and `3.11.2` release notes).
 - The GLM submodule has been removed, please use [Bambi](https://bambinos.github.io/bambi/) instead.
 - The `Distribution` keyword argument `testval` has been deprecated in favor of `initval`.
 - `pm.sample` now returns results as `InferenceData` instead of `MultiTrace` by default (see [#4744](https://github.com/pymc-devs/pymc3/pull/4744)).
 - `pm.sample_prior_predictive` no longer returns transformed variable values by default. Pass them by name in `var_names` if you want to obtain these draws (see [4769](https://github.com/pymc-devs/pymc3/pull/4769)).
-...
+- ...
 
 ### New Features
 - The `CAR` distribution has been added to allow for use of conditional autoregressions which often are used in spatial and network models.
@@ -20,12 +21,14 @@
 - Add `logcdf` method to Kumaraswamy distribution (see [#4706](https://github.com/pymc-devs/pymc3/pull/4706)).
 - The `OrderedMultinomial` distribution has been added for use on ordinal data which are _aggregated_ by trial, like multinomial observations, whereas `OrderedLogistic` only accepts ordinal data in a _disaggregated_ format, like categorical
   observations (see [#4773](https://github.com/pymc-devs/pymc3/pull/4773)).
+- ...
 
 ### Maintenance
 - Remove float128 dtype support (see [#4514](https://github.com/pymc-devs/pymc3/pull/4514)).
 - Logp method of `Uniform` and `DiscreteUniform` no longer depends on `pymc3.distributions.dist_math.bound` for proper evaluation (see [#4541](https://github.com/pymc-devs/pymc3/pull/4541)).
 - `Model.RV_dims` and `Model.coords` are now read-only properties. To modify the `coords` dictionary use `Model.add_coord`. Also `dims` or coordinate values that are `None` will be auto-completed (see [#4625](https://github.com/pymc-devs/pymc3/pull/4625)).
 - The length of `dims` in the model is now tracked symbolically through `Model.dim_lengths` (see [#4625](https://github.com/pymc-devs/pymc3/pull/4625)).
+- The `incomplete_beta` function in `pymc3.distributions.dist_math` was replaced by `aesara.tensor.betainc` (see [4736](https://github.com/pymc-devs/pymc3/pull/4736)).
 - ...
 
 ## PyMC3 3.11.2 (14 March 2021)

pymc3/distributions/continuous.py

@@ -59,7 +59,6 @@
     bound,
     clipped_beta_rvs,
     i0e,
-    incomplete_beta,
     log_i0,
     log_normal,
     logpow,
@@ -1246,23 +1245,19 @@ def logcdf(value, alpha, beta):
 
         Parameters
         ----------
-        value: numeric
-            Value(s) for which log CDF is calculated.
+        value: numeric or np.ndarray or aesara.tensor
+            Value(s) for which log CDF is calculated. If the log CDF for multiple
+            values are desired the values must be provided in a numpy array or Aesara tensor.
 
         Returns
         -------
         TensorVariable
         """
-        # incomplete_beta function can only handle scalar values (see #4342)
-        if np.ndim(value):
-            raise TypeError(
-                f"Beta.logcdf expects a scalar value but received a {np.ndim(value)}-dimensional object."
-            )
 
         return bound(
             at.switch(
                 at.lt(value, 1),
-                at.log(incomplete_beta(alpha, beta, value)),
+                at.log(at.betainc(alpha, beta, value)),
                 0,
             ),
             0 <= value,
@@ -1915,27 +1910,22 @@ def logcdf(value, nu, mu, sigma):
 
         Parameters
         ----------
-        value: numeric
-            Value(s) for which log CDF is calculated.
+        value: numeric or np.ndarray or aesara.tensor
+            Value(s) for which log CDF is calculated. If the log CDF for multiple
+            values are desired the values must be provided in a numpy array or Aesara tensor.
 
         Returns
         -------
         TensorVariable
         """
-        # incomplete_beta function can only handle scalar values (see #4342)
-        if np.ndim(value):
-            raise TypeError(
-                f"StudentT.logcdf expects a scalar value but received a {np.ndim(value)}-dimensional object."
-            )
-
         lam, sigma = get_tau_sigma(sigma=sigma)
 
         t = (value - mu) / sigma
         sqrt_t2_nu = at.sqrt(t ** 2 + nu)
         x = (t + sqrt_t2_nu) / (2.0 * sqrt_t2_nu)
 
         return bound(
-            at.log(incomplete_beta(nu / 2.0, nu / 2.0, x)),
+            at.log(at.betainc(nu / 2.0, nu / 2.0, x)),
             0 < nu,
             0 < sigma,
             0 < lam,
@@ -2380,13 +2370,8 @@ def logcdf(value, alpha, inv_beta):
         """
         beta = at.inv(inv_beta)
 
-        # Avoid C-assertion when the gammainc function is called with invalid values (#4340)
-        safe_alpha = at.switch(at.lt(alpha, 0), 0, alpha)
-        safe_beta = at.switch(at.lt(beta, 0), 0, beta)
-        safe_value = at.switch(at.lt(value, 0), 0, value)
-
         return bound(
-            at.log(at.gammainc(safe_alpha, safe_beta * safe_value)),
+            at.log(at.gammainc(alpha, beta * value)),
             0 <= value,
             0 < alpha,
             0 < beta,
@@ -2528,13 +2513,9 @@ def logcdf(value, alpha, beta):
         -------
         TensorVariable
         """
-        # Avoid C-assertion when the gammaincc function is called with invalid values (#4340)
-        safe_alpha = at.switch(at.lt(alpha, 0), 0, alpha)
-        safe_beta = at.switch(at.lt(beta, 0), 0, beta)
-        safe_value = at.switch(at.lt(value, 0), 0, value)
 
         return bound(
-            at.log(at.gammaincc(safe_alpha, safe_beta / safe_value)),
+            at.log(at.gammaincc(alpha, beta / value)),
             0 <= value,
             0 < alpha,
             0 < beta,

pymc3/distributions/discrete.py

@@ -35,7 +35,6 @@
     binomln,
     bound,
     factln,
-    incomplete_beta,
     log_diff_normal_cdf,
     logpow,
     normal_lccdf,
@@ -145,25 +144,20 @@ def logcdf(value, n, p):
 
         Parameters
         ----------
-        value: numeric
-            Value for which log CDF is calculated.
+        value: numeric or np.ndarray or aesara.tensor
+            Value(s) for which log CDF is calculated. If the log CDF for multiple
+            values are desired the values must be provided in a numpy array or Aesara tensor.
 
         Returns
         -------
         TensorVariable
         """
-        # incomplete_beta function can only handle scalar values (see #4342)
-        if np.ndim(value):
-            raise TypeError(
-                f"Binomial.logcdf expects a scalar value but received a {np.ndim(value)}-dimensional object."
-            )
-
         value = at.floor(value)
 
         return bound(
             at.switch(
                 at.lt(value, n),
-                at.log(incomplete_beta(n - value, value + 1, 1 - p)),
+                at.log(at.betainc(n - value, value + 1, 1 - p)),
                 0,
             ),
             0 <= value,
@@ -732,21 +726,16 @@ def logcdf(value, n, p):
 
         Parameters
         ----------
-        value: numeric
-            Value for which log CDF is calculated.
+        value: numeric or np.ndarray or aesara.tensor
+            Value(s) for which log CDF is calculated. If the log CDF for multiple
+            values are desired the values must be provided in a numpy array or Aesara tensor.
 
         Returns
         -------
         TensorVariable
         """
-        # incomplete_beta function can only handle scalar values (see #4342)
-        if np.ndim(value):
-            raise TypeError(
-                f"NegativeBinomial.logcdf expects a scalar value but received a {np.ndim(value)}-dimensional object."
-            )
-
         return bound(
-            at.log(incomplete_beta(n, at.floor(value) + 1, p)),
+            at.log(at.betainc(n, at.floor(value) + 1, p)),
             0 <= value,
             0 < n,
             0 <= p,
@@ -1461,20 +1450,15 @@ def logcdf(value, psi, n, p):
 
         Parameters
         ----------
-        value: numeric
-            Value for which log CDF is calculated.
+        value: numeric or np.ndarray or aesara.tensor
+            Value(s) for which log CDF is calculated. If the log CDF for multiple
+            values are desired the values must be provided in a numpy array or Aesara tensor.
 
         Returns
         -------
         TensorVariable
         """
 
-        # logcdf can only handle scalar values due to limitation in Binomial.logcdf
-        if np.ndim(value):
-            raise TypeError(
-                f"ZeroInflatedBinomial.logcdf expects a scalar value but received a {np.ndim(value)}-dimensional object."
-            )
-
         return bound(
             logaddexp(at.log1p(-psi), at.log(psi) + _logcdf(binomial, value, {}, n, p)),
             0 <= value,
@@ -1616,19 +1600,14 @@ def logcdf(value, psi, n, p):
 
         Parameters
         ----------
-        value: numeric
-            Value for which log CDF is calculated.
+        value: numeric or np.ndarray or aesara.tensor
+            Value(s) for which log CDF is calculated. If the log CDF for multiple
+            values are desired the values must be provided in a numpy array or Aesara tensor.
 
         Returns
         -------
         TensorVariable
         """
-        # logcdf can only handle scalar values due to limitation in NegativeBinomial.logcdf
-        if np.ndim(value):
-            raise TypeError(
-                f"ZeroInflatedNegativeBinomial.logcdf expects a scalar value but received a {np.ndim(value)}-dimensional object."
-            )
-
         return bound(
             logaddexp(at.log1p(-psi), at.log(psi) + _logcdf(nbinom, value, {}, n, p)),
             0 <= value,