diff --git a/src/pyrecest/distributions/nonperiodic/abstract_linear_distribution.py b/src/pyrecest/distributions/nonperiodic/abstract_linear_distribution.py
index 9b96b3037..7bd50d69e 100644
--- a/src/pyrecest/distributions/nonperiodic/abstract_linear_distribution.py
+++ b/src/pyrecest/distributions/nonperiodic/abstract_linear_distribution.py
@@ -9,6 +9,7 @@
 
 # pylint: disable=redefined-builtin,no-name-in-module,no-member
 # pylint: disable=no-name-in-module,no-member
+from pyrecest.backend import any as backend_any
 from pyrecest.backend import (
     array,
     atleast_1d,
@@ -316,6 +317,10 @@ def integrate_fun_over_domain(f, dim, left, right):
         right = AbstractLinearDistribution._normalize_static_integration_bound(
             right, dim, "right"
         )
+        if bool(backend_any(left > right)):
+            raise ValueError(
+                "left integration bound must not exceed right integration bound."
+            )
 
         def f_for_nquad(*args):
             # Avoid DeprecationWarning: Conversion of an array with ndim > 0 to a scalar is deprecated, and will error in future.
diff --git a/tests/distributions/test_abstract_linear_distribution.py b/tests/distributions/test_abstract_linear_distribution.py
index cc736f164..c1a09455d 100644
--- a/tests/distributions/test_abstract_linear_distribution.py
+++ b/tests/distributions/test_abstract_linear_distribution.py
@@ -63,6 +63,16 @@ def test_integrate_rejects_wrong_bound_shape(self):
         with self.assertRaisesRegex(ValueError, "right integration bound"):
             dist.integrate_numerically([0.0, 0.0], [1.0])
 
+    def test_integrate_rejects_reversed_bounds(self):
+        """Test that probability integration never accepts reversed intervals."""
+        dist_1d = GaussianDistribution(array([0.0]), array([[1.0]]))
+        with self.assertRaisesRegex(ValueError, "left integration bound"):
+            dist_1d.integrate_numerically(1.0, -1.0)
+
+        dist_2d = GaussianDistribution(array([0.0, 0.0]), diag(array([1.0, 1.0])))
+        with self.assertRaisesRegex(ValueError, "left integration bound"):
+            dist_2d.integrate_numerically([-1.0, 1.0], [1.0, -1.0])
+
     def test_integrate_fun_over_domain_1d_accepts_sequence_bounds(self):
         dist = GaussianDistribution(array([0.0]), array([[1.0]]))
 
@@ -91,6 +101,12 @@ def test_integrate_fun_over_domain_rejects_wrong_bound_shape(self):
                 lambda x: x[0], 2, [0.0], [1.0, 1.0]
             )
 
+    def test_integrate_fun_over_domain_rejects_reversed_bounds(self):
+        with self.assertRaisesRegex(ValueError, "left integration bound"):
+            AbstractLinearDistribution.integrate_fun_over_domain(
+                lambda x: x[0], 1, 1.0, -1.0
+            )
+
     def test_integrate_fun_over_domain(self):
         dist = GaussianDistribution(array([1.0, 2.0]), diag(array([1.0, 2.0])))