diff --git a/.github/workflows/pytest.yml b/.github/workflows/pytest.yml
index 1c18049d28..fad89b08f4 100644
--- a/.github/workflows/pytest.yml
+++ b/.github/workflows/pytest.yml
@@ -26,8 +26,6 @@ jobs:
         # 5th block: These tests PASS without a single XFAIL
         # 6th block: These have some XFAILs
           - |
-            --ignore=pymc3/tests/test_distribution_defaults.py
-            --ignore=pymc3/tests/test_distributions_random.py
             --ignore=pymc3/tests/test_distributions_timeseries.py
             --ignore=pymc3/tests/test_missing.py
             --ignore=pymc3/tests/test_mixture.py
@@ -53,7 +51,6 @@ jobs:
             --ignore=pymc3/tests/test_plots.py
             --ignore=pymc3/tests/test_special_functions.py
             --ignore=pymc3/tests/test_updates.py
-            --ignore=pymc3/tests/test_dist_math.py
             --ignore=pymc3/tests/test_examples.py
             --ignore=pymc3/tests/test_glm.py
             --ignore=pymc3/tests/test_gp.py
@@ -66,8 +63,11 @@ jobs:
             --ignore=pymc3/tests/test_posdef_sym.py
             --ignore=pymc3/tests/test_quadpotential.py
             --ignore=pymc3/tests/test_shape_handling.py
+            --ignore=pymc3/tests/test_distributions.py
+            --ignore=pymc3/tests/test_distributions_random.py
 
           - |
+            pymc3/tests/test_modelcontext.py
             pymc3/tests/test_dist_math.py
             pymc3/tests/test_minibatches.py
             pymc3/tests/test_pickling.py
@@ -76,7 +76,8 @@ jobs:
             pymc3/tests/test_updates.py
 
           - |
-            pymc3/tests/test_dist_math.py
+            pymc3/tests/test_distributions.py
+            pymc3/tests/test_distributions_random.py
             pymc3/tests/test_examples.py
             pymc3/tests/test_glm.py
             pymc3/tests/test_gp.py
diff --git a/pymc3/distributions/continuous.py b/pymc3/distributions/continuous.py
index 350358d71a..7377f1c4f3 100644
--- a/pymc3/distributions/continuous.py
+++ b/pymc3/distributions/continuous.py
@@ -23,7 +23,7 @@
 
 from aesara.assert_op import Assert
 from aesara.tensor.random.basic import (
-    beta,
+    BetaRV,
     cauchy,
     exponential,
     gamma,
@@ -42,6 +42,7 @@
     SplineWrapper,
     betaln,
     bound,
+    clipped_beta_rvs,
     gammaln,
     i0e,
     incomplete_beta,
@@ -786,18 +787,12 @@ def dist(cls, sigma=None, tau=None, sd=None, *args, **kwargs):
 
         tau, sigma = get_tau_sigma(tau=tau, sigma=sigma)
 
-        # sigma = sd = sigma = aet.as_tensor_variable(sigma)
-        # tau = tau = aet.as_tensor_variable(tau)
-
-        # mean = aet.sqrt(2 / (np.pi * tau))
-        # variance = (1.0 - 2 / np.pi) / tau
-
         assert_negative_support(tau, "tau", "HalfNormal")
         assert_negative_support(sigma, "sigma", "HalfNormal")
 
-        return super().dist([sigma, tau], **kwargs)
+        return super().dist([0.0, sigma], **kwargs)
 
-    def logp(value, sigma, tau):
+    def logp(value, loc, sigma):
         """
         Calculate log-probability of HalfNormal distribution at specified value.
 
@@ -811,14 +806,16 @@ def logp(value, sigma, tau):
         -------
         TensorVariable
         """
+        tau, sigma = get_tau_sigma(tau=None, sigma=sigma)
+
         return bound(
-            -0.5 * tau * value ** 2 + 0.5 * aet.log(tau * 2.0 / np.pi),
-            value >= 0,
+            -0.5 * tau * (value - loc) ** 2 + 0.5 * aet.log(tau * 2.0 / np.pi),
+            value >= loc,
             tau > 0,
             sigma > 0,
         )
 
-    def logcdf(value, sigma, tau):
+    def logcdf(value, loc, sigma):
         """
         Compute the log of the cumulative distribution function for HalfNormal distribution
         at the specified value.
@@ -833,10 +830,10 @@ def logcdf(value, sigma, tau):
         -------
         TensorVariable
         """
-        z = zvalue(value, mu=0, sigma=sigma)
+        z = zvalue(value, mu=loc, sigma=sigma)
         return bound(
             aet.log1p(-aet.erfc(z / aet.sqrt(2.0))),
-            0 <= value,
+            loc <= value,
             0 < sigma,
         )
 
@@ -1079,6 +1076,15 @@ def logcdf(self, value):
         )
 
 
+class BetaClippedRV(BetaRV):
+    @classmethod
+    def rng_fn(cls, rng, alpha, beta, size):
+        return clipped_beta_rvs(alpha, beta, size=size, random_state=rng)
+
+
+beta = BetaClippedRV()
+
+
 class Beta(UnitContinuous):
     r"""
     Beta log-likelihood.
@@ -1153,9 +1159,6 @@ def dist(cls, alpha=None, beta=None, mu=None, sigma=None, sd=None, *args, **kwar
         alpha = aet.as_tensor_variable(floatX(alpha))
         beta = aet.as_tensor_variable(floatX(beta))
 
-        # mean = alpha / (alpha + beta)
-        # variance = (alpha * beta) / ((alpha + beta) ** 2 * (alpha + beta + 1))
-
         assert_negative_support(alpha, "alpha", "Beta")
         assert_negative_support(beta, "beta", "Beta")
 
@@ -2259,14 +2262,10 @@ class HalfCauchy(PositiveContinuous):
     @classmethod
     def dist(cls, beta, *args, **kwargs):
         beta = aet.as_tensor_variable(floatX(beta))
-
-        # mode = aet.as_tensor_variable(0)
-        # median = beta
-
         assert_negative_support(beta, "beta", "HalfCauchy")
-        return super().dist([beta], **kwargs)
+        return super().dist([0.0, beta], **kwargs)
 
-    def logp(value, beta, alpha):
+    def logp(value, loc, beta):
         """
         Calculate log-probability of HalfCauchy distribution at specified value.
 
@@ -2281,12 +2280,12 @@ def logp(value, beta, alpha):
         TensorVariable
         """
         return bound(
-            aet.log(2) - aet.log(np.pi) - aet.log(beta) - aet.log1p((value / beta) ** 2),
-            value >= 0,
+            aet.log(2) - aet.log(np.pi) - aet.log(beta) - aet.log1p(((value - loc) / beta) ** 2),
+            value >= loc,
             beta > 0,
         )
 
-    def logcdf(value, beta, alpha):
+    def logcdf(value, loc, beta):
         """
         Compute the log of the cumulative distribution function for HalfCauchy distribution
         at the specified value.
@@ -2302,8 +2301,8 @@ def logcdf(value, beta, alpha):
         TensorVariable
         """
         return bound(
-            aet.log(2 * aet.arctan(value / beta) / np.pi),
-            0 <= value,
+            aet.log(2 * aet.arctan((value - loc) / beta) / np.pi),
+            loc <= value,
             0 < beta,
         )
 
diff --git a/pymc3/distributions/dist_math.py b/pymc3/distributions/dist_math.py
index 9bbe9fd30f..107dc95233 100644
--- a/pymc3/distributions/dist_math.py
+++ b/pymc3/distributions/dist_math.py
@@ -276,7 +276,7 @@ def MvNormalLogp():
     n, k = delta.shape
     n, k = f(n), f(k)
     chol_cov = cholesky(cov)
-    diag = aet.nlinalg.diag(chol_cov)
+    diag = aet.diag(chol_cov)
     ok = aet.all(diag > 0)
 
     chol_cov = aet.switch(ok, chol_cov, aet.fill(chol_cov, 1))
@@ -296,7 +296,7 @@ def dlogp(inputs, gradients):
         n, k = delta.shape
 
         chol_cov = cholesky(cov)
-        diag = aet.nlinalg.diag(chol_cov)
+        diag = aet.diag(chol_cov)
         ok = aet.all(diag > 0)
 
         chol_cov = aet.switch(ok, chol_cov, aet.fill(chol_cov, 1))
@@ -600,7 +600,7 @@ def incomplete_beta(a, b, value):
     )
 
 
-def clipped_beta_rvs(a, b, size=None, dtype="float64"):
+def clipped_beta_rvs(a, b, size=None, random_state=None, dtype="float64"):
     """Draw beta distributed random samples in the open :math:`(0, 1)` interval.
 
     The samples are generated with ``scipy.stats.beta.rvs``, but any value that
@@ -635,6 +635,6 @@ def clipped_beta_rvs(a, b, size=None, dtype="float64"):
         is shifted to ``np.nextafter(1, 0, dtype=dtype)``.
 
     """
-    out = scipy.stats.beta.rvs(a, b, size=size).astype(dtype)
+    out = scipy.stats.beta.rvs(a, b, size=size, random_state=random_state).astype(dtype)
     lower, upper = _beta_clip_values[dtype]
     return np.maximum(np.minimum(out, upper), lower)
diff --git a/pymc3/tests/test_dist_math.py b/pymc3/tests/test_dist_math.py
index c529c427f9..5ad725d6d6 100644
--- a/pymc3/tests/test_dist_math.py
+++ b/pymc3/tests/test_dist_math.py
@@ -19,6 +19,7 @@
 import numpy.testing as npt
 import pytest
 
+from aesara.tensor.random.basic import multinomial
 from scipy import interpolate, stats
 
 import pymc3 as pm
@@ -91,16 +92,13 @@ def test_alltrue_shape():
 
 
 class MultinomialA(Discrete):
-    def __init__(self, n, p, *args, **kwargs):
-        super().__init__(*args, **kwargs)
+    rv_op = multinomial
 
-        self.n = n
-        self.p = p
-
-    def logp(self, value):
-        n = self.n
-        p = self.p
+    @classmethod
+    def dist(cls, n, p, *args, **kwargs):
+        return super().dist([n, p], **kwargs)
 
+    def logp(value, n, p):
         return bound(
             factln(n) - factln(value).sum() + (value * aet.log(p)).sum(),
             value >= 0,
@@ -112,16 +110,13 @@ def logp(self, value):
 
 
 class MultinomialB(Discrete):
-    def __init__(self, n, p, *args, **kwargs):
-        super().__init__(*args, **kwargs)
-
-        self.n = n
-        self.p = p
+    rv_op = multinomial
 
-    def logp(self, value):
-        n = self.n
-        p = self.p
+    @classmethod
+    def dist(cls, n, p, *args, **kwargs):
+        return super().dist([n, p], **kwargs)
 
+    def logp(value, n, p):
         return bound(
             factln(n) - factln(value).sum() + (value * aet.log(p)).sum(),
             aet.all(value >= 0),
@@ -132,18 +127,17 @@ def logp(self, value):
         )
 
 
-@pytest.mark.xfail(reason="This test relies on the deprecated Distribution interface")
 def test_multinomial_bound():
 
     x = np.array([1, 5])
     n = x.sum()
 
     with pm.Model() as modelA:
-        p_a = pm.Dirichlet("p", floatX(np.ones(2)), shape=(2,))
+        p_a = pm.Dirichlet("p", floatX(np.ones(2)))
         MultinomialA("x", n, p_a, observed=x)
 
     with pm.Model() as modelB:
-        p_b = pm.Dirichlet("p", floatX(np.ones(2)), shape=(2,))
+        p_b = pm.Dirichlet("p", floatX(np.ones(2)))
         MultinomialB("x", n, p_b, observed=x)
 
     assert np.isclose(
@@ -151,7 +145,6 @@ def test_multinomial_bound():
     )
 
 
-@pytest.mark.xfail(reason="MvNormal not implemented")
 class TestMvNormalLogp:
     def test_logp(self):
         np.random.seed(42)
@@ -192,11 +185,10 @@ def func(chol_vec, delta):
         delta_val = floatX(np.random.randn(5, 2))
         verify_grad(func, [chol_vec_val, delta_val])
 
-    @pytest.mark.skip(reason="Fix in aesara not released yet: Theano#5908")
     @aesara.config.change_flags(compute_test_value="ignore")
     def test_hessian(self):
         chol_vec = aet.vector("chol_vec")
-        chol_vec.tag.test_value = np.array([0.1, 2, 3])
+        chol_vec.tag.test_value = floatX(np.array([0.1, 2, 3]))
         chol = aet.stack(
             [
                 aet.stack([aet.exp(0.1 * chol_vec[0]), 0]),
@@ -205,9 +197,10 @@ def test_hessian(self):
         )
         cov = aet.dot(chol, chol.T)
         delta = aet.matrix("delta")
-        delta.tag.test_value = np.ones((5, 2))
+        delta.tag.test_value = floatX(np.ones((5, 2)))
         logp = MvNormalLogp()(cov, delta)
         g_cov, g_delta = aet.grad(logp, [cov, delta])
+        # TODO: What's the test?  Something needs to be asserted.
         aet.grad(g_delta.sum() + g_cov.sum(), [delta, cov])
 
 
diff --git a/pymc3/tests/test_distribution_defaults.py b/pymc3/tests/test_distribution_defaults.py
deleted file mode 100644
index 4d0ecfe8b2..0000000000
--- a/pymc3/tests/test_distribution_defaults.py
+++ /dev/null
@@ -1,92 +0,0 @@
-#   Copyright 2020 The PyMC Developers
-#
-#   Licensed under the Apache License, Version 2.0 (the "License");
-#   you may not use this file except in compliance with the License.
-#   You may obtain a copy of the License at
-#
-#       http://www.apache.org/licenses/LICENSE-2.0
-#
-#   Unless required by applicable law or agreed to in writing, software
-#   distributed under the License is distributed on an "AS IS" BASIS,
-#   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
-#   See the License for the specific language governing permissions and
-#   limitations under the License.
-
-import numpy as np
-import pytest
-
-from pymc3.distributions import Categorical, Continuous, DiscreteUniform
-from pymc3.model import Model
-
-pytestmark = pytest.mark.xfail(reason="This test relies on the deprecated Distribution interface")
-
-
-class DistTest(Continuous):
-    def __init__(self, a, b, *args, **kwargs):
-        super().__init__(*args, **kwargs)
-        self.a = a
-        self.b = b
-
-    def logp(self, v):
-        return 0
-
-
-def test_default_nan_fail():
-    with Model(), pytest.raises(AttributeError):
-        DistTest("x", np.nan, 2, defaults=["a"])
-
-
-def test_default_empty_fail():
-    with Model(), pytest.raises(AttributeError):
-        DistTest("x", 1, 2, defaults=[])
-
-
-def test_default_testval():
-    with Model():
-        x = DistTest("x", 1, 2, testval=5, defaults=[])
-        assert x.tag.test_value == 5
-
-
-def test_default_testval_nan():
-    with Model():
-        x = DistTest("x", 1, 2, testval=np.nan, defaults=["a"])
-        np.testing.assert_almost_equal(x.tag.test_value, np.nan)
-
-
-def test_default_a():
-    with Model():
-        x = DistTest("x", 1, 2, defaults=["a"])
-        assert x.tag.test_value == 1
-
-
-def test_default_b():
-    with Model():
-        x = DistTest("x", np.nan, 2, defaults=["a", "b"])
-        assert x.tag.test_value == 2
-
-
-def test_default_c():
-    with Model():
-        y = DistTest("y", 7, 8, testval=94)
-        x = DistTest("x", y, 2, defaults=["a", "b"])
-        assert x.tag.test_value == 94
-
-
-def test_default_discrete_uniform():
-    with Model():
-        x = DiscreteUniform("x", lower=1, upper=2)
-        assert x.init_value == 1
-
-
-def test_discrete_uniform_negative():
-    model = Model()
-    with model:
-        x = DiscreteUniform("x", lower=-10, upper=0)
-    assert model.test_point["x"] == -5
-
-
-def test_categorical_mode():
-    model = Model()
-    with model:
-        x = Categorical("x", p=np.eye(4), shape=4)
-    assert np.allclose(model.test_point["x"], np.arange(4))
diff --git a/pymc3/tests/test_distributions_random.py b/pymc3/tests/test_distributions_random.py
index 95f379836a..19fdb6369b 100644
--- a/pymc3/tests/test_distributions_random.py
+++ b/pymc3/tests/test_distributions_random.py
@@ -28,6 +28,8 @@
 
 import pymc3 as pm
 
+from pymc3.aesaraf import floatX, intX
+from pymc3.distributions import change_rv_size
 from pymc3.distributions.dist_math import clipped_beta_rvs
 from pymc3.distributions.shape_utils import to_tuple
 from pymc3.exceptions import ShapeError
@@ -54,10 +56,6 @@
     product,
 )
 
-# XXX: This test module will need to be repurposed as tests for new
-# `RandomVariable`s and their `RandomVariable.perform` methods.
-pytestmark = pytest.mark.xfail(reason="This test relies on the deprecated Distribution interface")
-
 
 def pymc3_random(
     dist,
@@ -72,37 +70,64 @@ def pymc3_random(
 ):
     if model_args is None:
         model_args = {}
-    model = build_model(dist, valuedomain, paramdomains, extra_args)
+
+    model, param_vars = build_model(dist, valuedomain, paramdomains, extra_args)
+    model_dist = change_rv_size(model.named_vars["value"], size, expand=True)
+    pymc_rand = aesara.function([], model_dist)
+
     domains = paramdomains.copy()
     for pt in product(domains, n_samples=100):
         pt = pm.Point(pt, model=model)
         pt.update(model_args)
+
+        # Update the shared parameter variables in `param_vars`
+        for k, v in pt.items():
+            nv = param_vars.get(k, model.named_vars.get(k))
+            if nv.name in param_vars:
+                param_vars[nv.name].set_value(v)
+
         p = alpha
         # Allow KS test to fail (i.e., the samples be different)
         # a certain number of times. Crude, but necessary.
         f = fails
         while p <= alpha and f > 0:
-            s0 = model.named_vars["value"].random(size=size, point=pt)
-            s1 = ref_rand(size=size, **pt)
+            s0 = pymc_rand()
+            s1 = floatX(ref_rand(size=size, **pt))
             _, p = st.ks_2samp(np.atleast_1d(s0).flatten(), np.atleast_1d(s1).flatten())
             f -= 1
         assert p > alpha, str(pt)
 
 
 def pymc3_random_discrete(
-    dist, paramdomains, valuedomain=Domain([0]), ref_rand=None, size=100000, alpha=0.05, fails=20
+    dist,
+    paramdomains,
+    valuedomain=Domain([0]),
+    ref_rand=None,
+    size=100000,
+    alpha=0.05,
+    fails=20,
 ):
-    model = build_model(dist, valuedomain, paramdomains)
+    model, param_vars = build_model(dist, valuedomain, paramdomains)
+    model_dist = change_rv_size(model.named_vars["value"], size, expand=True)
+    pymc_rand = aesara.function([], model_dist)
+
     domains = paramdomains.copy()
     for pt in product(domains, n_samples=100):
         pt = pm.Point(pt, model=model)
         p = alpha
+
+        # Update the shared parameter variables in `param_vars`
+        for k, v in pt.items():
+            nv = param_vars.get(k, model.named_vars.get(k))
+            if nv.name in param_vars:
+                param_vars[nv.name].set_value(v)
+
         # Allow Chisq test to fail (i.e., the samples be different)
         # a certain number of times.
         f = fails
         while p <= alpha and f > 0:
-            o = model.named_vars["value"].random(size=size, point=pt)
-            e = ref_rand(size=size, **pt)
+            o = pymc_rand()
+            e = intX(ref_rand(size=size, **pt))
             o = np.atleast_1d(o).flatten()
             e = np.atleast_1d(e).flatten()
             observed = dict(zip(*np.unique(o, return_counts=True)))
@@ -149,7 +174,12 @@ def get_random_variable(self, shape, with_vector_params=False, name=None):
                         # in the test case parametrization "None" means "no specified (default)"
                         return self.distribution(name, transform=None, **params)
                     else:
-                        return self.distribution(name, shape=shape, transform=None, **params)
+                        ndim_supp = self.distribution.rv_op.ndim_supp
+                        if ndim_supp == 0:
+                            size = shape
+                        else:
+                            size = shape[:-ndim_supp]
+                        return self.distribution(name, size=size, transform=None, **params)
                 except TypeError:
                     if np.sum(np.atleast_1d(shape)) == 0:
                         pytest.skip("Timeseries must have positive shape")
@@ -158,16 +188,10 @@ def get_random_variable(self, shape, with_vector_params=False, name=None):
         @staticmethod
         def sample_random_variable(random_variable, size):
             """ Draws samples from a RandomVariable using its .random() method. """
-            try:
-                if size is None:
-                    return random_variable.random()
-                else:
-                    return random_variable.random(size=size)
-            except AttributeError:
-                if size is None:
-                    return random_variable.distribution.random()
-                else:
-                    return random_variable.distribution.random(size=size)
+            if size is None:
+                return random_variable.eval()
+            else:
+                return change_rv_size(random_variable, size, expand=True).eval()
 
         @pytest.mark.parametrize("size", [None, (), 1, (1,), 5, (4, 5)], ids=str)
         @pytest.mark.parametrize("shape", [None, ()], ids=str)
@@ -215,58 +239,63 @@ def test_vector_params(self, shape, size):
                 expected == actual
             ), f"Sample size {size} from {shape}-shaped RV had shape {actual}. Expected: {expected}"
 
-        @pytest.mark.parametrize("shape", [-2, 0, (0,), (2, 0), (5, 0, 3)])
-        def test_shape_error_on_zero_shape_rv(self, shape):
-            with pytest.raises(ValueError, match="not allowed"):
-                self.get_random_variable(shape)
-
 
+@pytest.mark.xfail(reason="This distribution has not been refactored for v4")
 class TestGaussianRandomWalk(BaseTestCases.BaseTestCase):
     distribution = pm.GaussianRandomWalk
     params = {"mu": 1.0, "sigma": 1.0}
     default_shape = (1,)
 
 
+@pytest.mark.skip(reason="This test is covered by Aesara")
 class TestNormal(BaseTestCases.BaseTestCase):
     distribution = pm.Normal
     params = {"mu": 0.0, "tau": 1.0}
 
 
+@pytest.mark.xfail(reason="This distribution has not been refactored for v4")
 class TestTruncatedNormal(BaseTestCases.BaseTestCase):
     distribution = pm.TruncatedNormal
     params = {"mu": 0.0, "tau": 1.0, "lower": -0.5, "upper": 0.5}
 
 
+@pytest.mark.xfail(reason="This distribution has not been refactored for v4")
 class TestTruncatedNormalLower(BaseTestCases.BaseTestCase):
     distribution = pm.TruncatedNormal
     params = {"mu": 0.0, "tau": 1.0, "lower": -0.5}
 
 
+@pytest.mark.xfail(reason="This distribution has not been refactored for v4")
 class TestTruncatedNormalUpper(BaseTestCases.BaseTestCase):
     distribution = pm.TruncatedNormal
     params = {"mu": 0.0, "tau": 1.0, "upper": 0.5}
 
 
+@pytest.mark.xfail(reason="This distribution has not been refactored for v4")
 class TestSkewNormal(BaseTestCases.BaseTestCase):
     distribution = pm.SkewNormal
     params = {"mu": 0.0, "sigma": 1.0, "alpha": 5.0}
 
 
+@pytest.mark.skip(reason="This test is covered by Aesara")
 class TestHalfNormal(BaseTestCases.BaseTestCase):
     distribution = pm.HalfNormal
     params = {"tau": 1.0}
 
 
+@pytest.mark.skip(reason="This test is covered by Aesara")
 class TestUniform(BaseTestCases.BaseTestCase):
     distribution = pm.Uniform
     params = {"lower": 0.0, "upper": 1.0}
 
 
+@pytest.mark.xfail(reason="This distribution has not been refactored for v4")
 class TestTriangular(BaseTestCases.BaseTestCase):
     distribution = pm.Triangular
     params = {"c": 0.5, "lower": 0.0, "upper": 1.0}
 
 
+@pytest.mark.xfail(reason="This distribution has not been refactored for v4")
 class TestWald(BaseTestCases.BaseTestCase):
     distribution = pm.Wald
     params = {"mu": 1.0, "lam": 1.0, "alpha": 0.0}
@@ -277,166 +306,199 @@ class TestBeta(BaseTestCases.BaseTestCase):
     params = {"alpha": 1.0, "beta": 1.0}
 
 
+@pytest.mark.xfail(reason="This distribution has not been refactored for v4")
 class TestKumaraswamy(BaseTestCases.BaseTestCase):
     distribution = pm.Kumaraswamy
     params = {"a": 1.0, "b": 1.0}
 
 
+@pytest.mark.skip(reason="This test is covered by Aesara")
 class TestExponential(BaseTestCases.BaseTestCase):
     distribution = pm.Exponential
     params = {"lam": 1.0}
 
 
+@pytest.mark.xfail(reason="This distribution has not been refactored for v4")
 class TestLaplace(BaseTestCases.BaseTestCase):
     distribution = pm.Laplace
     params = {"mu": 1.0, "b": 1.0}
 
 
+@pytest.mark.xfail(reason="This distribution has not been refactored for v4")
 class TestAsymmetricLaplace(BaseTestCases.BaseTestCase):
     distribution = pm.AsymmetricLaplace
     params = {"kappa": 1.0, "b": 1.0, "mu": 0.0}
 
 
+@pytest.mark.xfail(reason="This distribution has not been refactored for v4")
 class TestLognormal(BaseTestCases.BaseTestCase):
     distribution = pm.Lognormal
     params = {"mu": 1.0, "tau": 1.0}
 
 
+@pytest.mark.xfail(reason="This distribution has not been refactored for v4")
 class TestStudentT(BaseTestCases.BaseTestCase):
     distribution = pm.StudentT
     params = {"nu": 5.0, "mu": 0.0, "lam": 1.0}
 
 
+@pytest.mark.xfail(reason="This distribution has not been refactored for v4")
 class TestPareto(BaseTestCases.BaseTestCase):
     distribution = pm.Pareto
     params = {"alpha": 0.5, "m": 1.0}
 
 
+@pytest.mark.skip(reason="This test is covered by Aesara")
 class TestCauchy(BaseTestCases.BaseTestCase):
     distribution = pm.Cauchy
     params = {"alpha": 1.0, "beta": 1.0}
 
 
+@pytest.mark.skip(reason="This test is covered by Aesara")
 class TestHalfCauchy(BaseTestCases.BaseTestCase):
     distribution = pm.HalfCauchy
     params = {"beta": 1.0}
 
 
+@pytest.mark.skip(reason="This test is covered by Aesara")
 class TestGamma(BaseTestCases.BaseTestCase):
     distribution = pm.Gamma
     params = {"alpha": 1.0, "beta": 1.0}
 
 
+@pytest.mark.skip(reason="This test is covered by Aesara")
 class TestInverseGamma(BaseTestCases.BaseTestCase):
     distribution = pm.InverseGamma
     params = {"alpha": 0.5, "beta": 0.5}
 
 
+@pytest.mark.xfail(reason="This distribution has not been refactored for v4")
 class TestChiSquared(BaseTestCases.BaseTestCase):
     distribution = pm.ChiSquared
     params = {"nu": 2.0}
 
 
+@pytest.mark.xfail(reason="This distribution has not been refactored for v4")
 class TestWeibull(BaseTestCases.BaseTestCase):
     distribution = pm.Weibull
     params = {"alpha": 1.0, "beta": 1.0}
 
 
+@pytest.mark.xfail(reason="This distribution has not been refactored for v4")
 class TestExGaussian(BaseTestCases.BaseTestCase):
     distribution = pm.ExGaussian
     params = {"mu": 0.0, "sigma": 1.0, "nu": 1.0}
 
 
+@pytest.mark.xfail(reason="This distribution has not been refactored for v4")
 class TestVonMises(BaseTestCases.BaseTestCase):
     distribution = pm.VonMises
     params = {"mu": 0.0, "kappa": 1.0}
 
 
+@pytest.mark.xfail(reason="This distribution has not been refactored for v4")
 class TestGumbel(BaseTestCases.BaseTestCase):
     distribution = pm.Gumbel
     params = {"mu": 0.0, "beta": 1.0}
 
 
+@pytest.mark.xfail(reason="This distribution has not been refactored for v4")
 class TestLogistic(BaseTestCases.BaseTestCase):
     distribution = pm.Logistic
     params = {"mu": 0.0, "s": 1.0}
 
 
+@pytest.mark.xfail(reason="This distribution has not been refactored for v4")
 class TestLogitNormal(BaseTestCases.BaseTestCase):
     distribution = pm.LogitNormal
     params = {"mu": 0.0, "sigma": 1.0}
 
 
+@pytest.mark.skip(reason="This test is covered by Aesara")
 class TestBinomial(BaseTestCases.BaseTestCase):
     distribution = pm.Binomial
     params = {"n": 5, "p": 0.5}
 
 
+@pytest.mark.xfail(reason="This distribution has not been refactored for v4")
 class TestBetaBinomial(BaseTestCases.BaseTestCase):
     distribution = pm.BetaBinomial
     params = {"n": 5, "alpha": 1.0, "beta": 1.0}
 
 
+@pytest.mark.skip(reason="This test is covered by Aesara")
 class TestBernoulli(BaseTestCases.BaseTestCase):
     distribution = pm.Bernoulli
     params = {"p": 0.5}
 
 
+@pytest.mark.xfail(reason="This distribution has not been refactored for v4")
 class TestDiscreteWeibull(BaseTestCases.BaseTestCase):
     distribution = pm.DiscreteWeibull
     params = {"q": 0.25, "beta": 2.0}
 
 
+@pytest.mark.skip(reason="This test is covered by Aesara")
 class TestPoisson(BaseTestCases.BaseTestCase):
     distribution = pm.Poisson
     params = {"mu": 1.0}
 
 
+@pytest.mark.skip(reason="This test is covered by Aesara")
 class TestNegativeBinomial(BaseTestCases.BaseTestCase):
     distribution = pm.NegativeBinomial
     params = {"mu": 1.0, "alpha": 1.0}
 
 
+@pytest.mark.xfail(reason="This distribution has not been refactored for v4")
 class TestConstant(BaseTestCases.BaseTestCase):
     distribution = pm.Constant
     params = {"c": 3}
 
 
+@pytest.mark.xfail(reason="This distribution has not been refactored for v4")
 class TestZeroInflatedPoisson(BaseTestCases.BaseTestCase):
     distribution = pm.ZeroInflatedPoisson
     params = {"theta": 1.0, "psi": 0.3}
 
 
+@pytest.mark.xfail(reason="This distribution has not been refactored for v4")
 class TestZeroInflatedNegativeBinomial(BaseTestCases.BaseTestCase):
     distribution = pm.ZeroInflatedNegativeBinomial
     params = {"mu": 1.0, "alpha": 1.0, "psi": 0.3}
 
 
+@pytest.mark.xfail(reason="This distribution has not been refactored for v4")
 class TestZeroInflatedBinomial(BaseTestCases.BaseTestCase):
     distribution = pm.ZeroInflatedBinomial
     params = {"n": 10, "p": 0.6, "psi": 0.3}
 
 
+@pytest.mark.xfail(reason="This distribution has not been refactored for v4")
 class TestDiscreteUniform(BaseTestCases.BaseTestCase):
     distribution = pm.DiscreteUniform
     params = {"lower": 0.0, "upper": 10.0}
 
 
+@pytest.mark.xfail(reason="This distribution has not been refactored for v4")
 class TestGeometric(BaseTestCases.BaseTestCase):
     distribution = pm.Geometric
     params = {"p": 0.5}
 
 
+@pytest.mark.xfail(reason="This distribution has not been refactored for v4")
 class TestHyperGeometric(BaseTestCases.BaseTestCase):
     distribution = pm.HyperGeometric
     params = {"N": 50, "k": 25, "n": 10}
 
 
+@pytest.mark.xfail(reason="This distribution has not been refactored for v4")
 class TestMoyal(BaseTestCases.BaseTestCase):
     distribution = pm.Moyal
     params = {"mu": 0.0, "sigma": 1.0}
 
 
+@pytest.mark.skip(reason="This test is covered by Aesara")
 class TestCategorical(BaseTestCases.BaseTestCase):
     distribution = pm.Categorical
     params = {"p": np.ones(BaseTestCases.BaseTestCase.shape)}
@@ -456,6 +518,7 @@ def test_probability_vector_shape(self):
         assert pm.Categorical.dist(p=p).random(size=4).shape == (4, 3, 7)
 
 
+@pytest.mark.skip(reason="This test is covered by Aesara")
 class TestDirichlet(SeededTest):
     @pytest.mark.parametrize(
         "shape, size",
@@ -475,6 +538,7 @@ def test_dirichlet_random_shape(self, shape, size):
 
 
 class TestScalarParameterSamples(SeededTest):
+    @pytest.mark.xfail(reason="This distribution has not been refactored for v4")
     def test_bounded(self):
         # A bit crude...
         BoundedNormal = pm.Bound(pm.Normal, upper=0)
@@ -484,18 +548,21 @@ def ref_rand(size, tau):
 
         pymc3_random(BoundedNormal, {"tau": Rplus}, ref_rand=ref_rand)
 
+    @pytest.mark.skip(reason="This test is covered by Aesara")
     def test_uniform(self):
         def ref_rand(size, lower, upper):
             return st.uniform.rvs(size=size, loc=lower, scale=upper - lower)
 
         pymc3_random(pm.Uniform, {"lower": -Rplus, "upper": Rplus}, ref_rand=ref_rand)
 
+    @pytest.mark.skip(reason="This test is covered by Aesara")
     def test_normal(self):
         def ref_rand(size, mu, sigma):
             return st.norm.rvs(size=size, loc=mu, scale=sigma)
 
         pymc3_random(pm.Normal, {"mu": R, "sigma": Rplus}, ref_rand=ref_rand)
 
+    @pytest.mark.xfail(reason="This distribution has not been refactored for v4")
     def test_truncated_normal(self):
         def ref_rand(size, mu, sigma, lower, upper):
             return st.truncnorm.rvs(
@@ -508,6 +575,7 @@ def ref_rand(size, mu, sigma, lower, upper):
             ref_rand=ref_rand,
         )
 
+    @pytest.mark.xfail(reason="This distribution has not been refactored for v4")
     def test_truncated_normal_lower(self):
         def ref_rand(size, mu, sigma, lower):
             return st.truncnorm.rvs((lower - mu) / sigma, np.inf, size=size, loc=mu, scale=sigma)
@@ -516,6 +584,7 @@ def ref_rand(size, mu, sigma, lower):
             pm.TruncatedNormal, {"mu": R, "sigma": Rplusbig, "lower": -Rplusbig}, ref_rand=ref_rand
         )
 
+    @pytest.mark.xfail(reason="This distribution has not been refactored for v4")
     def test_truncated_normal_upper(self):
         def ref_rand(size, mu, sigma, upper):
             return st.truncnorm.rvs(-np.inf, (upper - mu) / sigma, size=size, loc=mu, scale=sigma)
@@ -524,18 +593,21 @@ def ref_rand(size, mu, sigma, upper):
             pm.TruncatedNormal, {"mu": R, "sigma": Rplusbig, "upper": Rplusbig}, ref_rand=ref_rand
         )
 
+    @pytest.mark.xfail(reason="This distribution has not been refactored for v4")
     def test_skew_normal(self):
         def ref_rand(size, alpha, mu, sigma):
             return st.skewnorm.rvs(size=size, a=alpha, loc=mu, scale=sigma)
 
         pymc3_random(pm.SkewNormal, {"mu": R, "sigma": Rplus, "alpha": R}, ref_rand=ref_rand)
 
+    @pytest.mark.skip(reason="This test is covered by Aesara")
     def test_half_normal(self):
         def ref_rand(size, tau):
             return st.halfnorm.rvs(size=size, loc=0, scale=tau ** -0.5)
 
         pymc3_random(pm.HalfNormal, {"tau": Rplus}, ref_rand=ref_rand)
 
+    @pytest.mark.xfail(reason="This distribution has not been refactored for v4")
     def test_wald(self):
         # Cannot do anything too exciting as scipy wald is a
         # location-scale model of the *standard* wald with mu=1 and lam=1
@@ -548,24 +620,28 @@ def ref_rand(size, mu, lam, alpha):
             ref_rand=ref_rand,
         )
 
+    @pytest.mark.skip(reason="This test is covered by Aesara")
     def test_beta(self):
         def ref_rand(size, alpha, beta):
             return clipped_beta_rvs(a=alpha, b=beta, size=size)
 
         pymc3_random(pm.Beta, {"alpha": Rplus, "beta": Rplus}, ref_rand=ref_rand)
 
+    @pytest.mark.skip(reason="This test is covered by Aesara")
     def test_exponential(self):
         def ref_rand(size, lam):
             return nr.exponential(scale=1.0 / lam, size=size)
 
         pymc3_random(pm.Exponential, {"lam": Rplus}, ref_rand=ref_rand)
 
+    @pytest.mark.xfail(reason="This distribution has not been refactored for v4")
     def test_laplace(self):
         def ref_rand(size, mu, b):
             return st.laplace.rvs(mu, b, size=size)
 
         pymc3_random(pm.Laplace, {"mu": R, "b": Rplus}, ref_rand=ref_rand)
 
+    @pytest.mark.xfail(reason="This distribution has not been refactored for v4")
     def test_laplace_asymmetric(self):
         def ref_rand(size, kappa, b, mu):
             u = np.random.uniform(size=size)
@@ -577,66 +653,77 @@ def ref_rand(size, kappa, b, mu):
 
         pymc3_random(pm.AsymmetricLaplace, {"b": Rplus, "kappa": Rplus, "mu": R}, ref_rand=ref_rand)
 
+    @pytest.mark.xfail(reason="This distribution has not been refactored for v4")
     def test_lognormal(self):
         def ref_rand(size, mu, tau):
             return np.exp(mu + (tau ** -0.5) * st.norm.rvs(loc=0.0, scale=1.0, size=size))
 
         pymc3_random(pm.Lognormal, {"mu": R, "tau": Rplusbig}, ref_rand=ref_rand)
 
+    @pytest.mark.xfail(reason="This distribution has not been refactored for v4")
     def test_student_t(self):
         def ref_rand(size, nu, mu, lam):
             return st.t.rvs(nu, mu, lam ** -0.5, size=size)
 
         pymc3_random(pm.StudentT, {"nu": Rplus, "mu": R, "lam": Rplus}, ref_rand=ref_rand)
 
+    @pytest.mark.skip(reason="This test is covered by Aesara")
     def test_cauchy(self):
         def ref_rand(size, alpha, beta):
             return st.cauchy.rvs(alpha, beta, size=size)
 
         pymc3_random(pm.Cauchy, {"alpha": R, "beta": Rplusbig}, ref_rand=ref_rand)
 
+    @pytest.mark.skip(reason="This test is covered by Aesara")
     def test_half_cauchy(self):
         def ref_rand(size, beta):
             return st.halfcauchy.rvs(scale=beta, size=size)
 
         pymc3_random(pm.HalfCauchy, {"beta": Rplusbig}, ref_rand=ref_rand)
 
+    @pytest.mark.skip(reason="This test is covered by Aesara")
     def test_gamma_alpha_beta(self):
         def ref_rand(size, alpha, beta):
             return st.gamma.rvs(alpha, scale=1.0 / beta, size=size)
 
         pymc3_random(pm.Gamma, {"alpha": Rplusbig, "beta": Rplusbig}, ref_rand=ref_rand)
 
+    @pytest.mark.skip(reason="This test is covered by Aesara")
     def test_gamma_mu_sigma(self):
         def ref_rand(size, mu, sigma):
             return st.gamma.rvs(mu ** 2 / sigma ** 2, scale=sigma ** 2 / mu, size=size)
 
         pymc3_random(pm.Gamma, {"mu": Rplusbig, "sigma": Rplusbig}, ref_rand=ref_rand)
 
+    @pytest.mark.skip(reason="This test is covered by Aesara")
     def test_inverse_gamma(self):
         def ref_rand(size, alpha, beta):
             return st.invgamma.rvs(a=alpha, scale=beta, size=size)
 
         pymc3_random(pm.InverseGamma, {"alpha": Rplus, "beta": Rplus}, ref_rand=ref_rand)
 
+    @pytest.mark.xfail(reason="This distribution has not been refactored for v4")
     def test_pareto(self):
         def ref_rand(size, alpha, m):
             return st.pareto.rvs(alpha, scale=m, size=size)
 
         pymc3_random(pm.Pareto, {"alpha": Rplusbig, "m": Rplusbig}, ref_rand=ref_rand)
 
+    @pytest.mark.xfail(reason="This distribution has not been refactored for v4")
     def test_ex_gaussian(self):
         def ref_rand(size, mu, sigma, nu):
             return nr.normal(mu, sigma, size=size) + nr.exponential(scale=nu, size=size)
 
         pymc3_random(pm.ExGaussian, {"mu": R, "sigma": Rplus, "nu": Rplus}, ref_rand=ref_rand)
 
+    @pytest.mark.xfail(reason="This distribution has not been refactored for v4")
     def test_vonmises(self):
         def ref_rand(size, mu, kappa):
             return st.vonmises.rvs(size=size, loc=mu, kappa=kappa)
 
         pymc3_random(pm.VonMises, {"mu": R, "kappa": Rplus}, ref_rand=ref_rand)
 
+    @pytest.mark.xfail(reason="This distribution has not been refactored for v4")
     def test_triangular(self):
         def ref_rand(size, lower, upper, c):
             scale = upper - lower
@@ -647,21 +734,25 @@ def ref_rand(size, lower, upper, c):
             pm.Triangular, {"lower": Runif, "upper": Runif + 3, "c": Runif + 1}, ref_rand=ref_rand
         )
 
+    @pytest.mark.xfail(reason="This distribution has not been refactored for v4")
     def test_flat(self):
         with pm.Model():
             f = pm.Flat("f")
             with pytest.raises(ValueError):
                 f.random(1)
 
+    @pytest.mark.xfail(reason="This distribution has not been refactored for v4")
     def test_half_flat(self):
         with pm.Model():
             f = pm.HalfFlat("f")
             with pytest.raises(ValueError):
                 f.random(1)
 
+    @pytest.mark.skip(reason="This test is covered by Aesara")
     def test_binomial(self):
         pymc3_random_discrete(pm.Binomial, {"n": Nat, "p": Unit}, ref_rand=st.binom.rvs)
 
+    @pytest.mark.xfail(reason="This distribution has not been refactored for v4")
     @pytest.mark.xfail(
         sys.platform.startswith("win"),
         reason="Known issue: https://github.com/pymc-devs/pymc3/pull/4269",
@@ -674,14 +765,17 @@ def test_beta_binomial(self):
     def _beta_bin(self, n, alpha, beta, size=None):
         return st.binom.rvs(n, st.beta.rvs(a=alpha, b=beta, size=size))
 
+    @pytest.mark.skip(reason="This test is covered by Aesara")
     def test_bernoulli(self):
         pymc3_random_discrete(
             pm.Bernoulli, {"p": Unit}, ref_rand=lambda size, p=None: st.bernoulli.rvs(p, size=size)
         )
 
+    @pytest.mark.skip(reason="This test is covered by Aesara")
     def test_poisson(self):
         pymc3_random_discrete(pm.Poisson, {"mu": Rplusbig}, size=500, ref_rand=st.poisson.rvs)
 
+    @pytest.mark.skip(reason="This test is covered by Aesara")
     def test_negative_binomial(self):
         def ref_rand(size, alpha, mu):
             return st.nbinom.rvs(alpha, alpha / (mu + alpha), size=size)
@@ -694,9 +788,11 @@ def ref_rand(size, alpha, mu):
             ref_rand=ref_rand,
         )
 
+    @pytest.mark.xfail(reason="This distribution has not been refactored for v4")
     def test_geometric(self):
         pymc3_random_discrete(pm.Geometric, {"p": Unit}, size=500, fails=50, ref_rand=nr.geometric)
 
+    @pytest.mark.xfail(reason="This distribution has not been refactored for v4")
     def test_hypergeometric(self):
         def ref_rand(size, N, k, n):
             return st.hypergeom.rvs(M=N, n=k, N=n, size=size)
@@ -713,6 +809,7 @@ def ref_rand(size, N, k, n):
             ref_rand=ref_rand,
         )
 
+    @pytest.mark.xfail(reason="This distribution has not been refactored for v4")
     def test_discrete_uniform(self):
         def ref_rand(size, lower, upper):
             return st.randint.rvs(lower, upper + 1, size=size)
@@ -721,6 +818,7 @@ def ref_rand(size, lower, upper):
             pm.DiscreteUniform, {"lower": -NatSmall, "upper": NatSmall}, ref_rand=ref_rand
         )
 
+    @pytest.mark.xfail(reason="This distribution has not been refactored for v4")
     def test_discrete_weibull(self):
         def ref_rand(size, q, beta):
             u = np.random.uniform(size=size)
@@ -731,6 +829,7 @@ def ref_rand(size, q, beta):
             pm.DiscreteWeibull, {"q": Unit, "beta": Rplusdunif}, ref_rand=ref_rand
         )
 
+    @pytest.mark.skip(reason="This test is covered by Aesara")
     @pytest.mark.parametrize("s", [2, 3, 4])
     def test_categorical_random(self, s):
         def ref_rand(size, p):
@@ -738,12 +837,14 @@ def ref_rand(size, p):
 
         pymc3_random_discrete(pm.Categorical, {"p": Simplex(s)}, ref_rand=ref_rand)
 
+    @pytest.mark.xfail(reason="This distribution has not been refactored for v4")
     def test_constant_dist(self):
         def ref_rand(size, c):
             return c * np.ones(size, dtype=int)
 
         pymc3_random_discrete(pm.Constant, {"c": I}, ref_rand=ref_rand)
 
+    @pytest.mark.skip(reason="This test is covered by Aesara")
     def test_mv_normal(self):
         def ref_rand(size, mu, cov):
             return st.multivariate_normal.rvs(mean=mu, cov=cov, size=size)
@@ -788,6 +889,7 @@ def ref_rand_uchol(size, mu, chol):
                 extra_args={"lower": False},
             )
 
+    @pytest.mark.xfail(reason="This distribution has not been refactored for v4")
     def test_matrix_normal(self):
         def ref_rand(size, mu, rowcov, colcov):
             return st.matrix_normal.rvs(mean=mu, rowcov=rowcov, colcov=colcov, size=size)
@@ -851,6 +953,7 @@ def ref_rand_uchol(size, mu, rowchol, colchol):
                     ref_rand=ref_rand_chol_transpose,
                 )
 
+    @pytest.mark.xfail(reason="This distribution has not been refactored for v4")
     def test_kronecker_normal(self):
         def ref_rand(size, mu, covs, sigma):
             cov = pm.math.kronecker(covs[0], covs[1]).eval()
@@ -912,6 +1015,7 @@ def ref_rand_evd(size, mu, evds, sigma):
                 model_args=evd_args,
             )
 
+    @pytest.mark.xfail(reason="This distribution has not been refactored for v4")
     def test_mv_t(self):
         def ref_rand(size, nu, Sigma, mu):
             normal = st.multivariate_normal.rvs(cov=Sigma, size=size)
@@ -927,6 +1031,7 @@ def ref_rand(size, nu, Sigma, mu):
                 ref_rand=ref_rand,
             )
 
+    @pytest.mark.skip(reason="This test is covered by Aesara")
     def test_dirichlet(self):
         def ref_rand(size, a):
             return st.dirichlet.rvs(a, size=size)
@@ -940,6 +1045,7 @@ def ref_rand(size, a):
                 ref_rand=ref_rand,
             )
 
+    @pytest.mark.xfail(reason="This distribution has not been refactored for v4")
     def test_dirichlet_multinomial(self):
         def ref_rand(size, a, n):
             k = a.shape[-1]
@@ -959,6 +1065,7 @@ def ref_rand(size, a, n):
                 ref_rand=ref_rand,
             )
 
+    @pytest.mark.xfail(reason="This distribution has not been refactored for v4")
     @pytest.mark.parametrize(
         "a, shape, n",
         [
@@ -990,6 +1097,7 @@ def test_dirichlet_multinomial_shape(self, a, shape, n):
         assert to_tuple(samp1.shape) == (1, *shape_)
         assert to_tuple(samp2.shape) == (2, *shape_)
 
+    @pytest.mark.xfail(reason="This distribution has not been refactored for v4")
     @pytest.mark.parametrize(
         "n, a, shape, expectation",
         [
@@ -1004,6 +1112,7 @@ def test_dirichlet_multinomial_dist_ShapeError(self, n, a, shape, expectation):
         with expectation:
             m.random()
 
+    @pytest.mark.skip(reason="This test is covered by Aesara")
     def test_multinomial(self):
         def ref_rand(size, p, n):
             return nr.multinomial(pvals=p, n=n, size=size)
@@ -1017,30 +1126,35 @@ def ref_rand(size, p, n):
                 ref_rand=ref_rand,
             )
 
+    @pytest.mark.xfail(reason="This distribution has not been refactored for v4")
     def test_gumbel(self):
         def ref_rand(size, mu, beta):
             return st.gumbel_r.rvs(loc=mu, scale=beta, size=size)
 
         pymc3_random(pm.Gumbel, {"mu": R, "beta": Rplus}, ref_rand=ref_rand)
 
+    @pytest.mark.xfail(reason="This distribution has not been refactored for v4")
     def test_logistic(self):
         def ref_rand(size, mu, s):
             return st.logistic.rvs(loc=mu, scale=s, size=size)
 
         pymc3_random(pm.Logistic, {"mu": R, "s": Rplus}, ref_rand=ref_rand)
 
+    @pytest.mark.xfail(reason="This distribution has not been refactored for v4")
     def test_logitnormal(self):
         def ref_rand(size, mu, sigma):
             return expit(st.norm.rvs(loc=mu, scale=sigma, size=size))
 
         pymc3_random(pm.LogitNormal, {"mu": R, "sigma": Rplus}, ref_rand=ref_rand)
 
+    @pytest.mark.xfail(reason="This distribution has not been refactored for v4")
     def test_moyal(self):
         def ref_rand(size, mu, sigma):
             return st.moyal.rvs(loc=mu, scale=sigma, size=size)
 
         pymc3_random(pm.Moyal, {"mu": R, "sigma": Rplus}, ref_rand=ref_rand)
 
+    @pytest.mark.xfail(reason="This distribution has not been refactored for v4")
     @pytest.mark.xfail(condition=(aesara.config.floatX == "float32"), reason="Fails on float32")
     def test_interpolated(self):
         for mu in R.vals:
@@ -1057,6 +1171,7 @@ def __init__(self, **kwargs):
 
                 pymc3_random(TestedInterpolated, {}, ref_rand=ref_rand)
 
+    @pytest.mark.xfail(reason="This distribution has not been refactored for v4")
     @pytest.mark.skip(
         "Wishart random sampling not implemented.\n"
         "See https://github.com/pymc-devs/pymc3/issues/538"
@@ -1072,6 +1187,7 @@ def test_wishart(self):
         #                           st.wishart(V, df=n, size=size))
         pass
 
+    @pytest.mark.xfail(reason="This distribution has not been refactored for v4")
     def test_lkj(self):
         for n in [2, 10, 50]:
             # pylint: disable=cell-var-from-loop
@@ -1093,6 +1209,7 @@ def __init__(self, **kwargs):
                 ref_rand=ref_rand,
             )
 
+    @pytest.mark.xfail(reason="This distribution has not been refactored for v4")
     def test_normalmixture(self):
         def ref_rand(size, w, mu, sigma):
             component = np.random.choice(w.size, size=size, p=w)
@@ -1122,6 +1239,7 @@ def ref_rand(size, w, mu, sigma):
         )
 
 
+@pytest.mark.xfail(reason="This distribution has not been refactored for v4")
 def test_mixture_random_shape():
     # test the shape broadcasting in mixture random
     y = np.concatenate([nr.poisson(5, size=10), nr.poisson(9, size=10)])
@@ -1159,7 +1277,7 @@ def test_mixture_random_shape():
     assert ppc["like3"].shape == (200, 20)
 
 
-@pytest.mark.xfail
+@pytest.mark.xfail(reason="This distribution has not been refactored for v4")
 def test_mixture_random_shape_fast():
     # test the shape broadcasting in mixture random
     y = np.concatenate([nr.poisson(5, size=10), nr.poisson(9, size=10)])
@@ -1190,6 +1308,7 @@ def test_mixture_random_shape_fast():
     assert rand3.shape == (100, 20)
 
 
+@pytest.mark.xfail(reason="This distribution has not been refactored for v4")
 class TestDensityDist:
     @pytest.mark.parametrize("shape", [(), (3,), (3, 2)], ids=str)
     def test_density_dist_with_random_sampleable(self, shape):
@@ -1299,6 +1418,7 @@ def test_density_dist_without_random_not_sampleable(self):
             pm.sample_posterior_predictive(trace, samples=samples, model=model, size=100)
 
 
+@pytest.mark.xfail(reason="This distribution has not been refactored for v4")
 class TestNestedRandom(SeededTest):
     def build_model(self, distribution, shape, nested_rvs_info):
         with pm.Model() as model:
@@ -1607,6 +1727,7 @@ def generate_shapes(include_params=False):
     return data
 
 
+@pytest.mark.skip(reason="This test is covered by Aesara")
 class TestMvNormal(SeededTest):
     @pytest.mark.parametrize(
         ["sample_shape", "dist_shape", "mu_shape", "param"],
@@ -1691,6 +1812,7 @@ def test_issue_3706(self):
         assert prior_pred["X"].shape == (1, N, 2)
 
 
+@pytest.mark.xfail(reason="This distribution has not been refactored for v4")
 def test_matrix_normal_random_with_random_variables():
     """
     This test checks for shape correctness when using MatrixNormal distribution
@@ -1714,6 +1836,7 @@ def test_matrix_normal_random_with_random_variables():
     assert prior["mu"].shape == (2, D, K)
 
 
+@pytest.mark.xfail(reason="This distribution has not been refactored for v4")
 class TestMvGaussianRandomWalk(SeededTest):
     @pytest.mark.parametrize(
         ["sample_shape", "dist_shape", "mu_shape", "param"],