diff --git a/RELEASE-NOTES.md b/RELEASE-NOTES.md
index 2cbff44fc3..f755ec145f 100644
--- a/RELEASE-NOTES.md
+++ b/RELEASE-NOTES.md
@@ -19,6 +19,7 @@ It also brings some dreadfully awaited fixes, so be sure to go through the chang
 - Option to set `check_bounds=False` when instantiating `pymc3.Model()`. This turns off bounds checks that ensure that input parameters of distributions are valid. For correctly specified models, this is unneccessary as all parameters get automatically transformed so that all values are valid. Turning this off should lead to faster sampling (see [#4377](https://github.com/pymc-devs/pymc3/pull/4377)).
 - `OrderedProbit` distribution added (see [#4232](https://github.com/pymc-devs/pymc3/pull/4232)).
 - `plot_posterior_predictive_glm` now works with `arviz.InferenceData` as well (see [#4234](https://github.com/pymc-devs/pymc3/pull/4234))
+- Add `logcdf` method to all univariate discrete distributions (see [#4387](https://github.com/pymc-devs/pymc3/pull/4387)).
 
 ### Maintenance
 - Fixed bug whereby partial traces returns after keyboard interrupt during parallel sampling had fewer draws than would've been available [#4318](https://github.com/pymc-devs/pymc3/pull/4318)
diff --git a/pymc3/distributions/discrete.py b/pymc3/distributions/discrete.py
index e639ba6684..d1dc858f9c 100644
--- a/pymc3/distributions/discrete.py
+++ b/pymc3/distributions/discrete.py
@@ -24,6 +24,7 @@
     binomln,
     bound,
     factln,
+    incomplete_beta,
     log_diff_normal_cdf,
     logpow,
     normal_lccdf,
@@ -32,7 +33,7 @@
 )
 from pymc3.distributions.distribution import Discrete, draw_values, generate_samples
 from pymc3.distributions.shape_utils import broadcast_distribution_samples
-from pymc3.math import log1pexp, logaddexp, logit, sigmoid, tround
+from pymc3.math import log1mexp, log1pexp, logaddexp, logit, logsumexp, sigmoid, tround
 from pymc3.theanof import floatX, intX, take_along_axis
 
 __all__ = [
@@ -148,6 +149,44 @@ def logp(self, value):
             p <= 1,
         )
 
+    def logcdf(self, value):
+        """
+        Compute the log of the cumulative distribution function for Binomial distribution
+        at the specified value.
+
+        Parameters
+        ----------
+        value: numeric
+            Value for which log CDF is calculated.
+
+        Returns
+        -------
+        TensorVariable
+        """
+        # incomplete_beta function can only handle scalar values (see #4342)
+        if np.ndim(value):
+            raise TypeError(
+                "Binomial.logcdf expects a scalar value but received a {}-dimensional object.".format(
+                    np.ndim(value)
+                )
+            )
+
+        n = self.n
+        p = self.p
+        value = tt.floor(value)
+
+        return bound(
+            tt.switch(
+                tt.lt(value, n),
+                tt.log(incomplete_beta(n - value, value + 1, 1 - p)),
+                0,
+            ),
+            0 <= value,
+            0 < n,
+            0 <= p,
+            p <= 1,
+        )
+
 
 class BetaBinomial(Discrete):
     R"""
@@ -271,16 +310,53 @@ def logp(self, value):
         """
         alpha = self.alpha
         beta = self.beta
+        n = self.n
         return bound(
-            binomln(self.n, value)
-            + betaln(value + alpha, self.n - value + beta)
-            - betaln(alpha, beta),
+            binomln(n, value) + betaln(value + alpha, n - value + beta) - betaln(alpha, beta),
             value >= 0,
-            value <= self.n,
+            value <= n,
             alpha > 0,
             beta > 0,
         )
 
+    def logcdf(self, value):
+        """
+        Compute the log of the cumulative distribution function for BetaBinomial distribution
+        at the specified value.
+
+        Parameters
+        ----------
+        value: numeric
+            Value for which log CDF is calculated.
+
+        Returns
+        -------
+        TensorVariable
+        """
+        # logcdf can only handle scalar values at the moment
+        if np.ndim(value):
+            raise TypeError(
+                "BetaBinomial.logcdf expects a scalar value but received a {}-dimensional object.".format(
+                    np.ndim(value)
+                )
+            )
+
+        alpha = self.alpha
+        beta = self.beta
+        n = self.n
+        safe_lower = tt.switch(tt.lt(value, 0), value, 0)
+
+        return bound(
+            tt.switch(
+                tt.lt(value, n),
+                logsumexp(self.logp(tt.arange(safe_lower, value + 1)), keepdims=False),
+                0,
+            ),
+            0 <= value,
+            0 < alpha,
+            0 < beta,
+        )
+
 
 class Bernoulli(Discrete):
     R"""Bernoulli log-likelihood
@@ -380,6 +456,34 @@ def logp(self, value):
                 tt.switch(value, tt.log(p), tt.log(1 - p)), value >= 0, value <= 1, p >= 0, p <= 1
             )
 
+    def logcdf(self, value):
+        """
+        Compute the log of the cumulative distribution function for Bernoulli distribution
+        at the specified value.
+
+        Parameters
+        ----------
+        value: numeric
+            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 theano tensor.
+
+        Returns
+        -------
+        TensorVariable
+        """
+        p = self.p
+
+        return bound(
+            tt.switch(
+                tt.lt(value, 1),
+                tt.log1p(-p),
+                0,
+            ),
+            0 <= value,
+            0 <= p,
+            p <= 1,
+        )
+
     def _distr_parameters_for_repr(self):
         return ["p"]
 
@@ -426,36 +530,11 @@ def DiscreteWeibull(q, b, x):
     def __init__(self, q, beta, *args, **kwargs):
         super().__init__(*args, defaults=("median",), **kwargs)
 
-        self.q = q = tt.as_tensor_variable(floatX(q))
-        self.beta = beta = tt.as_tensor_variable(floatX(beta))
+        self.q = tt.as_tensor_variable(floatX(q))
+        self.beta = tt.as_tensor_variable(floatX(beta))
 
         self.median = self._ppf(0.5)
 
-    def logp(self, value):
-        r"""
-        Calculate log-probability of DiscreteWeibull distribution at specified value.
-
-        Parameters
-        ----------
-        value: numeric
-            Value(s) for which log-probability is calculated. If the log probabilities for multiple
-            values are desired the values must be provided in a numpy array or theano tensor
-
-        Returns
-        -------
-        TensorVariable
-        """
-        q = self.q
-        beta = self.beta
-
-        return bound(
-            tt.log(tt.power(q, tt.power(value, beta)) - tt.power(q, tt.power(value + 1, beta))),
-            0 <= value,
-            0 < q,
-            q < 1,
-            0 < beta,
-        )
-
     def _ppf(self, p):
         r"""
         The percentile point function (the inverse of the cumulative
@@ -492,6 +571,56 @@ def random(self, point=None, size=None):
 
         return generate_samples(self._random, q, beta, dist_shape=self.shape, size=size)
 
+    def logp(self, value):
+        r"""
+        Calculate log-probability of DiscreteWeibull distribution at specified value.
+
+        Parameters
+        ----------
+        value: numeric
+            Value(s) for which log-probability is calculated. If the log probabilities for multiple
+            values are desired the values must be provided in a numpy array or theano tensor
+
+        Returns
+        -------
+        TensorVariable
+        """
+        q = self.q
+        beta = self.beta
+        return bound(
+            tt.log(tt.power(q, tt.power(value, beta)) - tt.power(q, tt.power(value + 1, beta))),
+            0 <= value,
+            0 < q,
+            q < 1,
+            0 < beta,
+        )
+
+    def logcdf(self, value):
+        """
+        Compute the log of the cumulative distribution function for Discrete Weibull distribution
+        at the specified value.
+
+        Parameters
+        ----------
+        value: numeric
+            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 theano tensor.
+
+        Returns
+        -------
+        TensorVariable
+        """
+        q = self.q
+        beta = self.beta
+
+        return bound(
+            tt.log1p(-tt.power(q, tt.power(value + 1, beta))),
+            0 <= value,
+            0 < q,
+            q < 1,
+            0 < beta,
+        )
+
 
 class Poisson(Discrete):
     R"""
@@ -581,6 +710,33 @@ def logp(self, value):
         # Return zero when mu and value are both zero
         return tt.switch(tt.eq(mu, 0) * tt.eq(value, 0), 0, log_prob)
 
+    def logcdf(self, value):
+        """
+        Compute the log of the cumulative distribution function for Poisson distribution
+        at the specified value.
+
+        Parameters
+        ----------
+        value: numeric
+            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 theano tensor.
+
+        Returns
+        -------
+        TensorVariable
+        """
+        mu = self.mu
+        value = tt.floor(value)
+        # To avoid gammaincc C-assertion when given invalid values (#4340)
+        safe_mu = tt.switch(tt.lt(mu, 0), 0, mu)
+        safe_value = tt.switch(tt.lt(value, 0), 0, value)
+
+        return bound(
+            tt.log(tt.gammaincc(safe_value + 1, safe_mu)),
+            0 <= value,
+            0 <= mu,
+        )
+
 
 class NegativeBinomial(Discrete):
     R"""
@@ -737,6 +893,40 @@ def logp(self, value):
         # Return Poisson when alpha gets very large.
         return tt.switch(tt.gt(alpha, 1e10), Poisson.dist(self.mu).logp(value), negbinom)
 
+    def logcdf(self, value):
+        """
+        Compute the log of the cumulative distribution function for NegativeBinomial distribution
+        at the specified value.
+
+        Parameters
+        ----------
+        value: numeric
+            Value for which log CDF is calculated.
+
+        Returns
+        -------
+        TensorVariable
+        """
+        # incomplete_beta function can only handle scalar values (see #4342)
+        if np.ndim(value):
+            raise TypeError(
+                "NegativeBinomial.logcdf expects a scalar value but received a {}-dimensional object.".format(
+                    np.ndim(value)
+                )
+            )
+
+        # TODO: avoid `p` recomputation if distribution was defined in terms of `p`
+        alpha = self.alpha
+        p = alpha / (self.mu + alpha)
+
+        return bound(
+            tt.log(incomplete_beta(alpha, tt.floor(value) + 1, p)),
+            0 <= value,
+            0 < alpha,
+            0 <= p,
+            p <= 1,
+        )
+
     def _distr_parameters_for_repr(self):
         return self._param_type
 
@@ -820,6 +1010,30 @@ def logp(self, value):
         p = self.p
         return bound(tt.log(p) + logpow(1 - p, value - 1), 0 <= p, p <= 1, value >= 1)
 
+    def logcdf(self, value):
+        """
+        Compute the log of the cumulative distribution function for Geometric distribution
+        at the specified value.
+
+        Parameters
+        ----------
+        value: numeric
+            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 theano tensor.
+
+        Returns
+        -------
+        TensorVariable
+        """
+        p = self.p
+
+        return bound(
+            log1mexp(-tt.log1p(-p) * value),
+            0 <= value,
+            0 <= p,
+            p <= 1,
+        )
+
 
 class HyperGeometric(Discrete):
     R"""
@@ -935,6 +1149,48 @@ def logp(self, value):
         upper = tt.switch(tt.lt(k, n), k, n)
         return bound(result, lower <= value, value <= upper)
 
+    def logcdf(self, value):
+        """
+        Compute the log of the cumulative distribution function for HyperGeometric distribution
+        at the specified value.
+
+        Parameters
+        ----------
+        value: numeric
+            Value for which log CDF is calculated.
+
+        Returns
+        -------
+        TensorVariable
+        """
+        # logcdf can only handle scalar values at the moment
+        if np.ndim(value):
+            raise TypeError(
+                "BetaBinomial.logcdf expects a scalar value but received a {}-dimensional object.".format(
+                    np.ndim(value)
+                )
+            )
+
+        # TODO: Use lower upper in locgdf for smarter logsumexp?
+        N = self.N
+        n = self.n
+        k = self.k
+        safe_lower = tt.switch(tt.lt(value, 0), value, 0)
+
+        return bound(
+            tt.switch(
+                tt.lt(value, n),
+                logsumexp(self.logp(tt.arange(safe_lower, value + 1)), keepdims=False),
+                0,
+            ),
+            0 <= value,
+            0 < N,
+            0 <= k,
+            0 <= n,
+            k <= N,
+            n <= N,
+        )
+
 
 class DiscreteUniform(Discrete):
     R"""
@@ -1025,6 +1281,30 @@ def logp(self, value):
         lower = self.lower
         return bound(-tt.log(upper - lower + 1), lower <= value, value <= upper)
 
+    def logcdf(self, value):
+        """
+        Compute the log of the cumulative distribution function for Discrete uniform distribution
+        at the specified value.
+
+        Parameters
+        ----------
+        value: numeric
+            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 theano tensor.
+
+        Returns
+        -------
+        TensorVariable
+        """
+        upper = self.upper
+        lower = self.lower
+
+        return bound(
+            tt.log(tt.minimum(tt.floor(value), upper) - lower + 1) - tt.log(upper - lower + 1),
+            lower <= value,
+            lower <= upper,
+        )
+
 
 class Categorical(Discrete):
     R"""
@@ -1263,7 +1543,7 @@ class ZeroInflatedPoisson(Discrete):
     def __init__(self, psi, theta, *args, **kwargs):
         super().__init__(*args, **kwargs)
         self.theta = theta = tt.as_tensor_variable(floatX(theta))
-        self.psi = psi = tt.as_tensor_variable(floatX(psi))
+        self.psi = tt.as_tensor_variable(floatX(psi))
         self.pois = Poisson.dist(theta)
         self.mode = self.pois.mode
 
@@ -1314,6 +1594,30 @@ def logp(self, value):
 
         return bound(logp_val, 0 <= value, 0 <= psi, psi <= 1, 0 <= theta)
 
+    def logcdf(self, value):
+        """
+        Compute the log of the cumulative distribution function for ZeroInflatedPoisson distribution
+        at the specified value.
+
+        Parameters
+        ----------
+        value: numeric
+            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 theano tensor.
+
+        Returns
+        -------
+        TensorVariable
+        """
+        psi = self.psi
+
+        return bound(
+            logaddexp(tt.log1p(-psi), tt.log(psi) + self.pois.logcdf(value)),
+            0 <= value,
+            0 <= psi,
+            psi <= 1,
+        )
+
 
 class ZeroInflatedBinomial(Discrete):
     R"""
@@ -1422,6 +1726,37 @@ def logp(self, value):
 
         return bound(logp_val, 0 <= value, value <= n, 0 <= psi, psi <= 1, 0 <= p, p <= 1)
 
+    def logcdf(self, value):
+        """
+        Compute the log of the cumulative distribution function for ZeroInflatedBinomial distribution
+        at the specified value.
+
+        Parameters
+        ----------
+        value: numeric
+            Value for which log CDF is calculated.
+
+        Returns
+        -------
+        TensorVariable
+        """
+        # logcdf can only handle scalar values due to limitation in Binomial.logcdf
+        if np.ndim(value):
+            raise TypeError(
+                "ZeroInflatedBinomial.logcdf expects a scalar value but received a {}-dimensional object.".format(
+                    np.ndim(value)
+                )
+            )
+
+        psi = self.psi
+
+        return bound(
+            logaddexp(tt.log1p(-psi), tt.log(psi) + self.bin.logcdf(value)),
+            0 <= value,
+            0 <= psi,
+            psi <= 1,
+        )
+
 
 class ZeroInflatedNegativeBinomial(Discrete):
     R"""
@@ -1561,6 +1896,36 @@ def logp(self, value):
 
         return bound(logp_val, 0 <= value, 0 <= psi, psi <= 1, mu > 0, alpha > 0)
 
+    def logcdf(self, value):
+        """
+        Compute the log of the cumulative distribution function for ZeroInflatedNegativeBinomial distribution
+        at the specified value.
+
+        Parameters
+        ----------
+        value: numeric
+            Value for which log CDF is calculated.
+
+        Returns
+        -------
+        TensorVariable
+        """
+        # logcdf can only handle scalar values due to limitation in NegativeBinomial.logcdf
+        if np.ndim(value):
+            raise TypeError(
+                "ZeroInflatedNegativeBinomial.logcdf expects a scalar value but received a {}-dimensional object.".format(
+                    np.ndim(value)
+                )
+            )
+        psi = self.psi
+
+        return bound(
+            logaddexp(tt.log1p(-psi), tt.log(psi) + self.nb.logcdf(value)),
+            0 <= value,
+            0 <= psi,
+            psi <= 1,
+        )
+
 
 class OrderedLogistic(Categorical):
     R"""
diff --git a/pymc3/tests/test_distributions.py b/pymc3/tests/test_distributions.py
index 16ed273488..9a21e00c15 100644
--- a/pymc3/tests/test_distributions.py
+++ b/pymc3/tests/test_distributions.py
@@ -91,7 +91,7 @@
     ZeroInflatedPoisson,
     continuous,
 )
-from pymc3.math import kronecker
+from pymc3.math import kronecker, logsumexp
 from pymc3.model import Deterministic, Model, Point
 from pymc3.tests.helpers import SeededTest, select_by_precision
 from pymc3.theanof import floatX
@@ -575,6 +575,28 @@ def check_logcdf(
                 err_msg=str(pt),
             )
 
+    def check_selfconsistency_discrete_logcdf(
+        self, distribution, domain, paramdomains, decimal=None, n_samples=100
+    ):
+        """
+        Check that logcdf of discrete distributions matches sum of logps up to value
+        """
+        domains = paramdomains.copy()
+        domains["value"] = domain
+        if decimal is None:
+            decimal = select_by_precision(float64=6, float32=3)
+        for pt in product(domains, n_samples=n_samples):
+            params = dict(pt)
+            value = params.pop("value")
+            values = np.arange(domain.lower, value + 1)
+            dist = distribution.dist(**params)
+            assert_almost_equal(
+                dist.logcdf(value).tag.test_value,
+                logsumexp(dist.logp(values), keepdims=False).tag.test_value,
+                decimal=decimal,
+                err_msg=str(pt),
+            )
+
     def check_int_to_1(self, model, value, domain, paramdomains):
         pdf = model.fastfn(exp(model.logpt))
         for pt in product(paramdomains, n_samples=10):
@@ -642,6 +664,17 @@ def test_discrete_unif(self):
             {"lower": -Rplusdunif, "upper": Rplusdunif},
             lambda value, lower, upper: sp.randint.logpmf(value, lower, upper + 1),
         )
+        self.check_logcdf(
+            DiscreteUniform,
+            Rdunif,
+            {"lower": -Rplusdunif, "upper": Rplusdunif},
+            lambda value, lower, upper: sp.randint.logcdf(value, lower, upper + 1),
+        )
+        self.check_selfconsistency_discrete_logcdf(
+            DiscreteUniform,
+            Rdunif,
+            {"lower": -Rplusdunif, "upper": Rplusdunif},
+        )
 
     def test_flat(self):
         self.pymc3_matches_scipy(Flat, Runif, {}, lambda value: 0)
@@ -801,32 +834,95 @@ def test_exponential(self):
 
     def test_geometric(self):
         self.pymc3_matches_scipy(
-            Geometric, Nat, {"p": Unit}, lambda value, p: np.log(sp.geom.pmf(value, p))
+            Geometric,
+            Nat,
+            {"p": Unit},
+            lambda value, p: np.log(sp.geom.pmf(value, p)),
+        )
+        self.check_logcdf(
+            Geometric,
+            Nat,
+            {"p": Unit},
+            lambda value, p: sp.geom.logcdf(value, p),
+        )
+        self.check_selfconsistency_discrete_logcdf(
+            Geometric,
+            Nat,
+            {"p": Unit},
         )
 
     def test_hypergeometric(self):
         def modified_scipy_hypergeom_logpmf(value, N, k, n):
+            # Convert nan to -np.inf
             original_res = sp.hypergeom.logpmf(value, N, k, n)
             return original_res if not np.isnan(original_res) else -np.inf
 
+        def modified_scipy_hypergeom_logcdf(value, N, k, n):
+            # Convert nan to -np.inf
+            original_res = sp.hypergeom.logcdf(value, N, k, n)
+
+            # Correct for scipy bug in logcdf method (see https://github.com/scipy/scipy/issues/13280)
+            if not np.isnan(original_res):
+                pmfs = sp.hypergeom.logpmf(np.arange(value + 1), N, k, n)
+                if np.all(np.isnan(pmfs)):
+                    original_res = np.nan
+
+            return original_res if not np.isnan(original_res) else -np.inf
+
         self.pymc3_matches_scipy(
             HyperGeometric,
             Nat,
             {"N": NatSmall, "k": NatSmall, "n": NatSmall},
-            lambda value, N, k, n: modified_scipy_hypergeom_logpmf(value, N, k, n),
+            modified_scipy_hypergeom_logpmf,
+        )
+        self.check_logcdf(
+            HyperGeometric,
+            Nat,
+            {"N": NatSmall, "k": NatSmall, "n": NatSmall},
+            modified_scipy_hypergeom_logcdf,
+        )
+        self.check_selfconsistency_discrete_logcdf(
+            HyperGeometric,
+            Nat,
+            {"N": NatSmall, "k": NatSmall, "n": NatSmall},
         )
 
     def test_negative_binomial(self):
-        def test_fun(value, mu, alpha):
+        def scipy_mu_alpha_logpmf(value, mu, alpha):
             return sp.nbinom.logpmf(value, alpha, 1 - mu / (mu + alpha))
 
-        self.pymc3_matches_scipy(NegativeBinomial, Nat, {"mu": Rplus, "alpha": Rplus}, test_fun)
+        def scipy_mu_alpha_logcdf(value, mu, alpha):
+            return sp.nbinom.logcdf(value, alpha, 1 - mu / (mu + alpha))
+
+        self.pymc3_matches_scipy(
+            NegativeBinomial,
+            Nat,
+            {"mu": Rplus, "alpha": Rplus},
+            scipy_mu_alpha_logpmf,
+        )
         self.pymc3_matches_scipy(
             NegativeBinomial,
             Nat,
             {"p": Unit, "n": Rplus},
             lambda value, p, n: sp.nbinom.logpmf(value, n, p),
         )
+        self.check_logcdf(
+            NegativeBinomial,
+            Nat,
+            {"mu": Rplus, "alpha": Rplus},
+            scipy_mu_alpha_logcdf,
+        )
+        self.check_logcdf(
+            NegativeBinomial,
+            Nat,
+            {"p": Unit, "n": Rplus},
+            lambda value, p, n: sp.nbinom.logcdf(value, n, p),
+        )
+        self.check_selfconsistency_discrete_logcdf(
+            NegativeBinomial,
+            Nat,
+            {"mu": Rplus, "alpha": Rplus},
+        )
 
     @pytest.mark.parametrize(
         "mu, p, alpha, n, expected",
@@ -1024,11 +1120,43 @@ def test_binomial(self):
             {"n": NatSmall, "p": Unit},
             lambda value, n, p: sp.binom.logpmf(value, n, p),
         )
+        self.check_logcdf(
+            Binomial,
+            Nat,
+            {"n": NatSmall, "p": Unit},
+            lambda value, n, p: sp.binom.logcdf(value, n, p),
+        )
+        self.check_selfconsistency_discrete_logcdf(
+            Binomial,
+            Nat,
+            {"n": NatSmall, "p": Unit},
+        )
 
     # Too lazy to propagate decimal parameter through the whole chain of deps
     @pytest.mark.xfail(condition=(theano.config.floatX == "float32"), reason="Fails on float32")
     def test_beta_binomial(self):
-        self.checkd(BetaBinomial, Nat, {"alpha": Rplus, "beta": Rplus, "n": NatSmall})
+        self.checkd(
+            BetaBinomial,
+            Nat,
+            {"alpha": Rplus, "beta": Rplus, "n": NatSmall},
+        )
+        self.pymc3_matches_scipy(
+            BetaBinomial,
+            Nat,
+            {"alpha": Rplus, "beta": Rplus, "n": NatSmall},
+            lambda value, alpha, beta, n: sp.betabinom.logpmf(value, a=alpha, b=beta, n=n),
+        )
+        self.check_logcdf(
+            BetaBinomial,
+            Nat,
+            {"alpha": Rplus, "beta": Rplus, "n": NatSmall},
+            lambda value, alpha, beta, n: sp.betabinom.logcdf(value, a=alpha, b=beta, n=n),
+        )
+        self.check_selfconsistency_discrete_logcdf(
+            BetaBinomial,
+            Nat,
+            {"alpha": Rplus, "beta": Rplus, "n": NatSmall},
+        )
 
     def test_bernoulli(self):
         self.pymc3_matches_scipy(
@@ -1038,7 +1166,27 @@ def test_bernoulli(self):
             lambda value, logit_p: sp.bernoulli.logpmf(value, scipy.special.expit(logit_p)),
         )
         self.pymc3_matches_scipy(
-            Bernoulli, Bool, {"p": Unit}, lambda value, p: sp.bernoulli.logpmf(value, p)
+            Bernoulli,
+            Bool,
+            {"p": Unit},
+            lambda value, p: sp.bernoulli.logpmf(value, p),
+        )
+        self.check_logcdf(
+            Bernoulli,
+            Bool,
+            {"p": Unit},
+            lambda value, p: sp.bernoulli.logcdf(value, p),
+        )
+        self.check_logcdf(
+            Bernoulli,
+            Bool,
+            {"logit_p": R},
+            lambda value, logit_p: sp.bernoulli.logcdf(value, scipy.special.expit(logit_p)),
+        )
+        self.check_selfconsistency_discrete_logcdf(
+            Bernoulli,
+            Bool,
+            {"p": Unit},
         )
 
     def test_discrete_weibull(self):
@@ -1048,10 +1196,29 @@ def test_discrete_weibull(self):
             {"q": Unit, "beta": Rplusdunif},
             discrete_weibull_logpmf,
         )
+        self.check_selfconsistency_discrete_logcdf(
+            DiscreteWeibull,
+            Nat,
+            {"q": Unit, "beta": Rplusdunif},
+        )
 
     def test_poisson(self):
         self.pymc3_matches_scipy(
-            Poisson, Nat, {"mu": Rplus}, lambda value, mu: sp.poisson.logpmf(value, mu)
+            Poisson,
+            Nat,
+            {"mu": Rplus},
+            lambda value, mu: sp.poisson.logpmf(value, mu),
+        )
+        self.check_logcdf(
+            Poisson,
+            Nat,
+            {"mu": Rplus},
+            lambda value, mu: sp.poisson.logcdf(value, mu),
+        )
+        self.check_selfconsistency_discrete_logcdf(
+            Poisson,
+            Nat,
+            {"mu": Rplus},
         )
 
     def test_bound_poisson(self):
@@ -1073,7 +1240,16 @@ def test_constantdist(self):
     # Too lazy to propagate decimal parameter through the whole chain of deps
     @pytest.mark.xfail(condition=(theano.config.floatX == "float32"), reason="Fails on float32")
     def test_zeroinflatedpoisson(self):
-        self.checkd(ZeroInflatedPoisson, Nat, {"theta": Rplus, "psi": Unit})
+        self.checkd(
+            ZeroInflatedPoisson,
+            Nat,
+            {"theta": Rplus, "psi": Unit},
+        )
+        self.check_selfconsistency_discrete_logcdf(
+            ZeroInflatedPoisson,
+            Nat,
+            {"theta": Rplus, "psi": Unit},
+        )
 
     # Too lazy to propagate decimal parameter through the whole chain of deps
     @pytest.mark.xfail(condition=(theano.config.floatX == "float32"), reason="Fails on float32")
@@ -1083,11 +1259,25 @@ def test_zeroinflatednegativebinomial(self):
             Nat,
             {"mu": Rplusbig, "alpha": Rplusbig, "psi": Unit},
         )
+        self.check_selfconsistency_discrete_logcdf(
+            ZeroInflatedNegativeBinomial,
+            Nat,
+            {"mu": Rplusbig, "alpha": Rplusbig, "psi": Unit},
+        )
 
     # Too lazy to propagate decimal parameter through the whole chain of deps
     @pytest.mark.xfail(condition=(theano.config.floatX == "float32"), reason="Fails on float32")
     def test_zeroinflatedbinomial(self):
-        self.checkd(ZeroInflatedBinomial, Nat, {"n": NatSmall, "p": Unit, "psi": Unit})
+        self.checkd(
+            ZeroInflatedBinomial,
+            Nat,
+            {"n": NatSmall, "p": Unit, "psi": Unit},
+        )
+        self.check_selfconsistency_discrete_logcdf(
+            ZeroInflatedBinomial,
+            Nat,
+            {"n": NatSmall, "p": Unit, "psi": Unit},
+        )
 
     @pytest.mark.parametrize("n", [1, 2, 3])
     def test_mvnormal(self, n):