diff --git a/pymc/tests/test_distributions.py b/pymc/tests/test_distributions.py
index d29a929e97..2a9a4af219 100644
--- a/pymc/tests/test_distributions.py
+++ b/pymc/tests/test_distributions.py
@@ -35,6 +35,8 @@ def polyagamma_cdf(*args, **kwargs):
         raise RuntimeError("polyagamma package is not installed!")
 
 
+from contextlib import ExitStack as does_not_raise
+
 import pytest
 import scipy.stats
 import scipy.stats.distributions as sp
@@ -155,6 +157,14 @@ def __init__(self, vals, dtype=None, edges=None, shape=None):
         if edges is None:
             edges = array(vals[0]), array(vals[-1])
             vals = vals[1:-1]
+
+        if not vals:
+            raise ValueError(
+                f"Domain has no values left after removing edges: {edges}.\n"
+                "You can duplicate the edge values or explicitly specify the edges with the edge keyword.\n"
+                f"For example: `Domain([{edges[0]}, {edges[0]}, {edges[1]}, {edges[1]}])`"
+            )
+
         if shape is None:
             shape = avals[0].shape
 
@@ -192,6 +202,22 @@ def __neg__(self):
         return Domain([-v for v in self.vals], self.dtype, (-self.lower, -self.upper), self.shape)
 
 
+@pytest.mark.parametrize(
+    "values, edges, expectation",
+    [
+        ([], None, pytest.raises(IndexError)),
+        ([], (0, 0), pytest.raises(ValueError)),
+        ([0], None, pytest.raises(ValueError)),
+        ([0], (0, 0), does_not_raise()),
+        ([-1, 1], None, pytest.raises(ValueError)),
+        ([-1, 0, 1], None, does_not_raise()),
+    ],
+)
+def test_domain(values, edges, expectation):
+    with expectation:
+        Domain(values, edges=edges)
+
+
 def product(domains, n_samples=-1):
     """Get an iterator over a product of domains.
 
@@ -2423,7 +2449,7 @@ def test_categorical_valid_p(self):
     def test_categorical(self, n):
         self.check_logp(
             Categorical,
-            Domain(range(n), "int64"),
+            Domain(range(n), dtype="int64", edges=(None, None)),
             {"p": Simplex(n)},
             lambda value, p: categorical_logpdf(value, p),
         )
@@ -2432,7 +2458,7 @@ def test_categorical(self, n):
     def test_orderedlogistic(self, n):
         self.check_logp(
             OrderedLogistic,
-            Domain(range(n), "int64"),
+            Domain(range(n), dtype="int64", edges=(None, None)),
             {"eta": R, "cutpoints": Vector(R, n - 1)},
             lambda value, eta, cutpoints: orderedlogistic_logpdf(value, eta, cutpoints),
         )
@@ -2441,7 +2467,7 @@ def test_orderedlogistic(self, n):
     def test_orderedprobit(self, n):
         self.check_logp(
             OrderedProbit,
-            Domain(range(n), "int64"),
+            Domain(range(n), dtype="int64", edges=(None, None)),
             {"eta": Runif, "cutpoints": UnitSortedVector(n - 1)},
             lambda value, eta, cutpoints: orderedprobit_logpdf(value, eta, cutpoints),
         )
diff --git a/pymc/tests/test_distributions_random.py b/pymc/tests/test_distributions_random.py
index 8b7c189e6c..2797be839e 100644
--- a/pymc/tests/test_distributions_random.py
+++ b/pymc/tests/test_distributions_random.py
@@ -64,13 +64,16 @@ def pymc_random(
     dist,
     paramdomains,
     ref_rand,
-    valuedomain=Domain([0]),
+    valuedomain=None,
     size=10000,
     alpha=0.05,
     fails=10,
     extra_args=None,
     model_args=None,
 ):
+    if valuedomain is None:
+        valuedomain = Domain([0], edges=(None, None))
+
     if model_args is None:
         model_args = {}
 
@@ -104,12 +107,15 @@ def pymc_random(
 def pymc_random_discrete(
     dist,
     paramdomains,
-    valuedomain=Domain([0]),
+    valuedomain=None,
     ref_rand=None,
     size=100000,
     alpha=0.05,
     fails=20,
 ):
+    if valuedomain is None:
+        valuedomain = Domain([0], edges=(None, None))
+
     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)