diff --git a/docs/source/contributing/implementing_distribution.md b/docs/source/contributing/implementing_distribution.md
index 4d3ce8b4a6..38ab96d762 100644
--- a/docs/source/contributing/implementing_distribution.md
+++ b/docs/source/contributing/implementing_distribution.md
@@ -350,9 +350,11 @@ These methods will perform a grid evaluation on the combinations of `domain` and
 There are a couple of details worth keeping in mind:
 
 1. By default, the first and last values (edges) of the `Domain` are not compared (they are used for other things). If it is important to test the edge of the `Domain`, the edge values can be repeated. This is done by the `Bool`: `Bool = Domain([0, 0, 1, 1], "int64")`
-3. There are some default domains (such as `R` and `Rplus`) that you can use for testing your new distribution, but it's also perfectly fine to create your own domains inside the test function if there is a good reason for it (e.g., when the default values lead too many extreme unlikely combinations that are not very informative about the correctness of the implementation).
+3. There are some default domains (such as `R` and `Rplus`) that you can use for testing your new distribution, but it's also perfectly fine to create your own domains inside the test function if there is a good reason for it (e.g., when the default values lead too many extreme unlikely combinations that are not very informative about the correctness of the implementation). For example, there is an `Rplusbig` domain defined which is suited to scale parameters, where numbers close to zero are atypical.
 4. By default, a random subset of 100 `param` x `paramdomain` combinations is tested, to keep the test runtime under control. When testing your shiny new distribution, you can temporarily set `n_samples=-1` to force all combinations to be tested. This is important to avoid your `PR` leading to surprising failures in future runs whenever some bad combinations of parameters are randomly tested.
-5. On GitHub some tests run twice, under the `pytensor.config.floatX` flags of `"float64"` and `"float32"`. However, the reference Python functions will run in a pure "float64" environment, which means the reference and the PyMC results can diverge quite a lot (e.g., underflowing to `-np.inf` for extreme parameters). You should therefore make sure you test locally in both regimes. A quick and dirty way of doing this is to temporarily add `pytensor.config.floatX = "float32"` at the very top of file, immediately after `import pytensor`. Remember to set `n_samples=-1` as well to test all combinations. The test output will show what exact parameter values lead to a failure. If you are confident that your implementation is correct, you may opt to tweak the decimal precision with `select_by_precision`, or adjust the tested `Domain` values. In extreme cases, you can mark the test with a conditional `xfail` (if only one of the sub-methods is failing, they should be separated, so that the `xfail` is as narrow as possible):
+5. On GitHub some tests run twice, under the `pytensor.config.floatX` flags of `"float64"` and `"float32"`. However, the reference Python functions will run in a pure "float64" environment, which means the reference and the PyMC results can diverge quite a lot (e.g., underflowing to `-np.inf` for extreme parameters). You should therefore make sure you test locally in both regimes. A quick and dirty way of doing this is to temporarily add `pytensor.config.floatX = "float32"` at the very top of file, immediately after `import pytensor`. Remember to set `n_samples=-1` as well to test all combinations. The test output will show what exact parameter values lead to a failure. The tolerance (see next note) for `"float64"` or `"float32"` can be set with `select_by_precision`, or adjust the tested `Domain` values. In extreme cases, you can mark the test with a conditional `xfail` (if only one of the sub-methods is failing, they should be separated, so that the `xfail` is as narrow as possible):
+6. In the event of tests failing, and if you are confident that your implementation is correct, you may opt to tweak the tolerances against which your `logp` or `logc` is being tested. For this, you can pass an `atol` and/or `rtol` parmeter to `logp` or `logc`, which set the acceptable absolute and relative tolerances respectively, following the relevant {doc}`numpy:reference/generated/testing.assert_allclose <NumPy documentation>` for the `assert_allclose` which is utilized.
+
 
 ```python
 
diff --git a/pymc/distributions/continuous.py b/pymc/distributions/continuous.py
index dae67808d7..c902a4f571 100644
--- a/pymc/distributions/continuous.py
+++ b/pymc/distributions/continuous.py
@@ -12,7 +12,7 @@
 #   See the License for the specific language governing permissions and
 #   limitations under the License.
 
-# Contains code from Aeppl, Copyright (c) 2021-2022, PyTensor Developers.
+# Contains code from Aeppl, Copyright (c) 2021-2022, Aesara Developers.
 
 # coding: utf-8
 """
diff --git a/pymc/tests/distributions/test_continuous.py b/pymc/tests/distributions/test_continuous.py
index ba01ca0b62..5103655505 100644
--- a/pymc/tests/distributions/test_continuous.py
+++ b/pymc/tests/distributions/test_continuous.py
@@ -235,14 +235,14 @@ def test_polyagamma(self):
             Rplus,
             {"h": Rplus, "z": R},
             lambda value, h, z: polyagamma_pdf(value, h, z, return_log=True),
-            decimal=select_by_precision(float64=6, float32=-1),
+            rtol=select_by_precision(float64=1e-6, float32=1e-1),
         )
         check_logcdf(
             pm.PolyaGamma,
             Rplus,
             {"h": Rplus, "z": R},
             lambda value, h, z: polyagamma_cdf(value, h, z, return_log=True),
-            decimal=select_by_precision(float64=6, float32=-1),
+            rtol=select_by_precision(float64=1e-6, float32=1e-1),
         )
 
     def test_flat(self):
@@ -269,14 +269,14 @@ def test_normal(self):
             R,
             {"mu": R, "sigma": Rplus},
             lambda value, mu, sigma: st.norm.logpdf(value, mu, sigma),
-            decimal=select_by_precision(float64=6, float32=1),
+            rtol=select_by_precision(float64=1e-6, float32=1e-3),
         )
         check_logcdf(
             pm.Normal,
             R,
             {"mu": R, "sigma": Rplus},
             lambda value, mu, sigma: st.norm.logcdf(value, mu, sigma),
-            decimal=select_by_precision(float64=6, float32=1),
+            rtol=select_by_precision(float64=1e-6, float32=1e-3),
         )
 
     def test_half_normal(self):
@@ -285,7 +285,7 @@ def test_half_normal(self):
             Rplus,
             {"sigma": Rplus},
             lambda value, sigma: st.halfnorm.logpdf(value, scale=sigma),
-            decimal=select_by_precision(float64=6, float32=-1),
+            rtol=select_by_precision(float64=1e-6, float32=1e-1),
         )
         check_logcdf(
             pm.HalfNormal,
@@ -320,7 +320,7 @@ def test_wald_logp(self):
             Rplus,
             {"mu": Rplus, "alpha": Rplus},
             lambda value, mu, alpha: st.invgauss.logpdf(value, mu=mu, loc=alpha),
-            decimal=select_by_precision(float64=6, float32=1),
+            rtol=select_by_precision(float64=1e-6, float32=1e-2),
         )
 
     def test_wald_logcdf(self):
@@ -443,7 +443,7 @@ def test_laplace_asymmetric(self):
             R,
             {"b": Rplus, "kappa": Rplus, "mu": R},
             laplace_asymmetric_logpdf,
-            decimal=select_by_precision(float64=6, float32=2),
+            rtol=select_by_precision(float64=1e-6, float32=1e0),
         )
 
     def test_lognormal(self):
@@ -596,7 +596,7 @@ def test_fun(value, mu, sigma):
             Rplus,
             {"mu": Rplus, "sigma": Rplus},
             test_fun,
-            decimal=select_by_precision(float64=4, float32=3),
+            rtol=select_by_precision(float64=1e-6, float32=1e-5),
         )
 
     def test_pareto(self):
@@ -653,7 +653,7 @@ def test_skew_normal(self):
             R,
             {"mu": R, "sigma": Rplusbig, "alpha": R},
             lambda value, alpha, mu, sigma: st.skewnorm.logpdf(value, alpha, mu, sigma),
-            decimal=select_by_precision(float64=5, float32=3),
+            rtol=select_by_precision(float64=1e-7, float32=1e-5),
         )
 
     @pytest.mark.parametrize(
@@ -720,14 +720,14 @@ def test_logistic(self):
             R,
             {"mu": R, "s": Rplus},
             lambda value, mu, s: st.logistic.logpdf(value, mu, s),
-            decimal=select_by_precision(float64=6, float32=1),
+            rtol=select_by_precision(float64=1e-6, float32=1e-1),
         )
         check_logcdf(
             pm.Logistic,
             R,
             {"mu": R, "s": Rplus},
             lambda value, mu, s: st.logistic.logcdf(value, mu, s),
-            decimal=select_by_precision(float64=6, float32=1),
+            rtol=select_by_precision(float64=1e-6, float32=1e-1),
         )
 
     def test_logitnormal(self):
@@ -738,7 +738,7 @@ def test_logitnormal(self):
             lambda value, mu, sigma: (
                 st.norm.logpdf(sp.logit(value), mu, sigma) - (np.log(value) + np.log1p(-value))
             ),
-            decimal=select_by_precision(float64=6, float32=1),
+            rtol=select_by_precision(float64=1e-6, float32=1e-1),
         )
 
     @pytest.mark.skipif(
@@ -841,7 +841,7 @@ def scipy_logp(value, mu, sigma, lower, upper):
             R,
             {"mu": R, "sigma": Rplusbig, "lower": -Rplusbig, "upper": Rplusbig},
             scipy_logp,
-            decimal=select_by_precision(float64=6, float32=1),
+            rtol=select_by_precision(float64=1e-6, float32=1e-1),
             skip_paramdomain_outside_edge_test=True,
         )
 
@@ -850,7 +850,7 @@ def scipy_logp(value, mu, sigma, lower, upper):
             R,
             {"mu": R, "sigma": Rplusbig, "upper": Rplusbig},
             ft.partial(scipy_logp, lower=-np.inf),
-            decimal=select_by_precision(float64=6, float32=1),
+            rtol=select_by_precision(float64=1e-6, float32=1e-1),
             skip_paramdomain_outside_edge_test=True,
         )
 
@@ -859,7 +859,7 @@ def scipy_logp(value, mu, sigma, lower, upper):
             R,
             {"mu": R, "sigma": Rplusbig, "lower": -Rplusbig},
             ft.partial(scipy_logp, upper=np.inf),
-            decimal=select_by_precision(float64=6, float32=1),
+            rtol=select_by_precision(float64=1e-6, float32=1e-1),
             skip_paramdomain_outside_edge_test=True,
         )
 
diff --git a/pymc/tests/distributions/test_logprob.py b/pymc/tests/distributions/test_logprob.py
index 0133b4383a..0c0c13d2ca 100644
--- a/pymc/tests/distributions/test_logprob.py
+++ b/pymc/tests/distributions/test_logprob.py
@@ -224,7 +224,7 @@ def test_joint_logp_subtensor():
     # The compiled graph should not contain any `RandomVariables`
     assert_no_rvs(logp_vals_fn.maker.fgraph.outputs[0])
 
-    decimals = select_by_precision(float64=6, float32=4)
+    atol = select_by_precision(float64=1e-8, float32=1e-5)
 
     for i in range(10):
         bern_sp = sp.bernoulli(p)
@@ -238,7 +238,7 @@ def test_joint_logp_subtensor():
 
         logp_vals = logp_vals_fn(A_idx_value, I_value)
 
-        np.testing.assert_almost_equal(logp_vals, exp_obs_logps, decimal=decimals)
+        np.testing.assert_allclose(logp_vals, exp_obs_logps, atol=atol)
 
 
 def test_logp_helper():
diff --git a/pymc/tests/distributions/test_multivariate.py b/pymc/tests/distributions/test_multivariate.py
index f5b3fc04e3..ec62515219 100644
--- a/pymc/tests/distributions/test_multivariate.py
+++ b/pymc/tests/distributions/test_multivariate.py
@@ -272,7 +272,7 @@ def test_mvnormal(self, n):
             RealMatrix(5, n),
             {"mu": Vector(R, n), "chol": PdMatrixChol(n)},
             normal_logpdf_chol,
-            decimal=select_by_precision(float64=6, float32=-1),
+            rtol=select_by_precision(float64=1e-8, float32=1e-3),
             extra_args={"size": 5},
         )
         check_logp(
@@ -280,14 +280,14 @@ def test_mvnormal(self, n):
             Vector(R, n),
             {"mu": Vector(R, n), "chol": PdMatrixChol(n)},
             normal_logpdf_chol,
-            decimal=select_by_precision(float64=6, float32=0),
+            rtol=select_by_precision(float64=1e-8, float32=1e-1),
         )
         check_logp(
             pm.MvNormal,
             Vector(R, n),
             {"mu": Vector(R, n), "chol": PdMatrixCholUpper(n)},
             normal_logpdf_chol_upper,
-            decimal=select_by_precision(float64=6, float32=0),
+            rtol=select_by_precision(float64=1e-8, float32=1e-1),
             extra_args={"lower": False},
         )
 
@@ -338,7 +338,7 @@ def test_matrixnormal(self, n):
                 "colcov": PdMatrix(n) * mat_scale,
             },
             matrix_normal_logpdf_cov,
-            decimal=select_by_precision(float64=5, float32=3),
+            rtol=select_by_precision(float64=1e-8, float32=1e-6),
         )
         check_logp(
             pm.MatrixNormal,
@@ -349,7 +349,7 @@ def test_matrixnormal(self, n):
                 "colcov": PdMatrix(n) * mat_scale,
             },
             matrix_normal_logpdf_cov,
-            decimal=select_by_precision(float64=5, float32=3),
+            rtol=select_by_precision(float64=1e-8, float32=1e-6),
         )
         check_logp(
             pm.MatrixNormal,
@@ -360,7 +360,7 @@ def test_matrixnormal(self, n):
                 "colchol": PdMatrixChol(n) * mat_scale,
             },
             matrix_normal_logpdf_chol,
-            decimal=select_by_precision(float64=5, float32=3),
+            rtol=select_by_precision(float64=1e-8, float32=1e-6),
         )
         check_logp(
             pm.MatrixNormal,
@@ -371,7 +371,7 @@ def test_matrixnormal(self, n):
                 "colchol": PdMatrixChol(3) * mat_scale,
             },
             matrix_normal_logpdf_chol,
-            decimal=select_by_precision(float64=5, float32=3),
+            rtol=select_by_precision(float64=5, float32=3),
         )
 
     @pytest.mark.parametrize("n", [2, 3])
@@ -787,10 +787,8 @@ def test_car_logp(self, sparse, size):
         delta_logp = scipy_logp - car_logp
 
         # Check to make sure all the delta values are identical.
-        tol = 1e-08
-        if pytensor.config.floatX == "float32":
-            tol = 1e-5
-        assert np.allclose(delta_logp - delta_logp[0], 0.0, atol=tol)
+        atol = select_by_precision(float64=1e-8, float32=1e-5)
+        assert np.allclose(delta_logp - delta_logp[0], 0.0, atol=atol)
 
 
 @pytest.mark.parametrize(
diff --git a/pymc/tests/distributions/test_timeseries.py b/pymc/tests/distributions/test_timeseries.py
index 587adf1c2c..57cad14c57 100644
--- a/pymc/tests/distributions/test_timeseries.py
+++ b/pymc/tests/distributions/test_timeseries.py
@@ -754,8 +754,8 @@ def test_logp(self):
             z = Normal("z", mu=0, sigma=vol, shape=data.shape)
         garch_like = t.compile_logp(y)({"y": data})
         reg_like = t.compile_logp(z)({"z": data})
-        decimal = select_by_precision(float64=7, float32=4)
-        np.testing.assert_allclose(garch_like, reg_like, 10 ** (-decimal))
+        rtol = select_by_precision(float64=1e-9, float32=1e-6)
+        np.testing.assert_allclose(garch_like, reg_like, rtol)
 
     @pytest.mark.parametrize(
         "batched_param",
diff --git a/pymc/tests/distributions/util.py b/pymc/tests/distributions/util.py
index 6af649adb7..0fa3ceb692 100644
--- a/pymc/tests/distributions/util.py
+++ b/pymc/tests/distributions/util.py
@@ -232,7 +232,8 @@ def check_logp(
     domain,
     paramdomains,
     scipy_logp,
-    decimal=None,
+    rtol=None,
+    atol=None,
     n_samples=100,
     extra_args=None,
     scipy_args=None,
@@ -255,9 +256,14 @@ def check_logp(
         Supported domains of distribution parameters
     scipy_logp : Scipy logpmf/logpdf method
         Scipy logp method of equivalent pymc_dist distribution
-    decimal : Int
-        Level of precision with which pymc_dist and scipy logp are compared.
-        Defaults to 6 for float64 and 3 for float32
+    atol : Float
+        Absolute tolerance with which pymc_dist and scipy logp are comparing
+        the difference between actual and desired to `atol + rtol * abs(desired)`.
+        Defaults to 0
+    rtol : Float
+        Absolute tolerance with which pymc_dist and scipy logp are comparing
+        the difference between actual and desired to `atol + rtol * abs(desired)`.
+        Defaults to 1e-8 for float64 and 1e-5 for float32
     n_samples : Int
         Upper limit on the number of valid domain and value combinations that
         are compared between pymc and scipy methods. If n_samples is below the
@@ -269,8 +275,11 @@ def check_logp(
     scipy_args : Dictionary with extra arguments needed to call scipy logp method
         Usually the same as extra_args
     """
-    if decimal is None:
-        decimal = select_by_precision(float64=6, float32=3)
+    if rtol is None:
+        rtol = select_by_precision()
+
+    if atol is None:
+        atol = select_by_precision(float64=1e-10, float32=1e-8)
 
     if extra_args is None:
         extra_args = {}
@@ -317,10 +326,11 @@ def _model_input_dict(model, param_vars, pt):
         pt_d = _model_input_dict(model, param_vars, pt)
         pt_logp = pm.Point(pt_d, model=model)
         pt_ref = pm.Point(pt, filter_model_vars=False, model=model)
-        npt.assert_almost_equal(
+        npt.assert_allclose(
             logp_pymc(pt_logp),
             logp_reference(pt_ref),
-            decimal=decimal,
+            rtol=rtol,
+            atol=atol,
             err_msg=str(pt),
         )
 
@@ -383,7 +393,8 @@ def check_logcdf(
     domain,
     paramdomains,
     scipy_logcdf,
-    decimal=None,
+    rtol=None,
+    atol=None,
     n_samples=100,
     skip_paramdomain_inside_edge_test=False,
     skip_paramdomain_outside_edge_test=False,
@@ -413,9 +424,14 @@ def check_logcdf(
         Supported domains of distribution parameters
     scipy_logcdf : Scipy logcdf method
         Scipy logcdf method of equivalent pymc_dist distribution
-    decimal : Int
-        Level of precision with which pymc_dist and scipy_logcdf are compared.
-        Defaults to 6 for float64 and 3 for float32
+    atol : Float
+        Absolute tolerance with which pymc_dist and scipy logp are comparing
+        the difference between actual and desired to `atol + rtol * abs(desired)`.
+        Defaults to 0
+    rtol : Float
+        Absolute tolerance with which pymc_dist and scipy logp are comparing
+        the difference between actual and desired to `atol + rtol * abs(desired)`.
+        Defaults to 1e-8 for float64 and 1e-5 for float32
     n_samples : Int
         Upper limit on the number of valid domain and value combinations that
         are compared between pymc and scipy methods. If n_samples is below the
@@ -443,8 +459,11 @@ def check_logcdf(
         value = model.rvs_to_values[rv]
         pymc_logcdf = model.compile_fn(logcdf(rv, value))
 
-        if decimal is None:
-            decimal = select_by_precision(float64=6, float32=3)
+        if rtol is None:
+            rtol = select_by_precision()
+
+        if atol is None:
+            atol = select_by_precision(float64=1e-10, float32=1e-8)
 
         for pt in product(domains, n_samples=n_samples):
             params = dict(pt)
@@ -457,10 +476,11 @@ def check_logcdf(
             pymc_eval = pymc_logcdf({"value": value})
 
             params["value"] = value  # for displaying in err_msg
-            npt.assert_almost_equal(
+            npt.assert_allclose(
                 pymc_eval,
                 scipy_eval,
-                decimal=decimal,
+                rtol=rtol,
+                atol=atol,
                 err_msg=str(params),
             )
 
@@ -529,7 +549,8 @@ def check_selfconsistency_discrete_logcdf(
     distribution,
     domain,
     paramdomains,
-    decimal=None,
+    rtol=None,
+    atol=None,
     n_samples=100,
 ):
     """
@@ -537,8 +558,12 @@ def check_selfconsistency_discrete_logcdf(
     """
     domains = paramdomains.copy()
     domains["value"] = domain
-    if decimal is None:
-        decimal = select_by_precision(float64=6, float32=3)
+
+    if rtol is None:
+        rtol = select_by_precision(float64=1e-8, float32=1e-5)
+
+    if atol is None:
+        atol = select_by_precision(float64=1e-10, float32=1e-8)
 
     model, param_vars = build_model(distribution, domain, paramdomains)
     rv = model["value"]
@@ -556,10 +581,11 @@ def check_selfconsistency_discrete_logcdf(
             param_vars[param_name].set_value(param_value)
 
         with pytensor.config.change_flags(mode=Mode("py")):
-            npt.assert_almost_equal(
+            npt.assert_allclose(
                 dist_logcdf({"value": value}),
                 sp.logsumexp([dist_logp({"value": value}) for value in values]),
-                decimal=decimal,
+                rtol=rtol,
+                atol=atol,
                 err_msg=str(pt),
             )
 
diff --git a/pymc/tests/helpers.py b/pymc/tests/helpers.py
index ab0d7c7a4b..0b79241abb 100644
--- a/pymc/tests/helpers.py
+++ b/pymc/tests/helpers.py
@@ -119,10 +119,10 @@ def match_value(self, k, dv, v):
         return result
 
 
-def select_by_precision(float64, float32):
-    """Helper function to choose reasonable decimal cutoffs for different floatX modes."""
-    decimal = float64 if pytensor.config.floatX == "float64" else float32
-    return decimal
+def select_by_precision(float64=1e-6, float32=1e-3):
+    """Helper function to choose a precision tolerance for different floatX modes."""
+    tol = float64 if pytensor.config.floatX == "float64" else float32
+    return tol
 
 
 @contextlib.contextmanager