diff --git a/src/pyrecest/distributions/cart_prod/abstract_hypercylindrical_distribution.py b/src/pyrecest/distributions/cart_prod/abstract_hypercylindrical_distribution.py
index 8cbf58573..16b6a0909 100644
--- a/src/pyrecest/distributions/cart_prod/abstract_hypercylindrical_distribution.py
+++ b/src/pyrecest/distributions/cart_prod/abstract_hypercylindrical_distribution.py
@@ -29,6 +29,7 @@
     ndim,
     ones,
     pi,
+    reshape,
     sin,
     sqrt,
     tile,
@@ -63,6 +64,45 @@ def integrate_numerically(self, integration_boundaries=None):
         if integration_boundaries is None:
             integration_boundaries = self.get_reasonable_integration_boundaries()
 
+        try:
+            n_boundaries = len(integration_boundaries)
+        except TypeError as exc:
+            raise ValueError(
+                "integration_boundaries must contain one [lower, upper] interval "
+                "per input dimension"
+            ) from exc
+        if n_boundaries != self.input_dim:
+            raise ValueError(
+                "integration_boundaries must contain one [lower, upper] interval "
+                "per input dimension"
+            )
+        for boundary in integration_boundaries:
+            if callable(boundary):
+                continue
+            try:
+                interval = asarray(boundary)
+            except (TypeError, ValueError) as exc:
+                raise ValueError(
+                    "each integration boundary must be a [lower, upper] interval"
+                ) from exc
+            if interval.ndim != 1 or interval.shape[0] != 2:
+                raise ValueError(
+                    "each integration boundary must be a [lower, upper] interval"
+                )
+            try:
+                contains_nan = bool(any(isnan(interval)))
+                reversed_interval = bool(interval[0] > interval[1])
+            except (TypeError, ValueError) as exc:
+                raise ValueError(
+                    "integration boundaries must be real numeric intervals"
+                ) from exc
+            if contains_nan:
+                raise ValueError("integration boundaries must not contain NaN")
+            if reversed_interval:
+                raise ValueError(
+                    "integration boundary lower bound must not exceed upper bound"
+                )
+
         def f(*args):
             return self.pdf(array(args))
 
@@ -209,6 +249,7 @@ def f_cond_unnorm(xs, input_lin=input_lin):
             elif xs.ndim == 1:
                 if self.bound_dim != xs.shape[0]:
                     raise ValueError("Input should be of size (bound_dim,).")
+                xs = reshape(xs, (1, self.bound_dim))
                 n_inputs = 1
             else:
                 n_inputs = xs.shape[0]
@@ -261,6 +302,7 @@ def f_cond_unnorm(xs, input_periodic=input_periodic):
             elif xs.ndim == 1:
                 if self.lin_dim != xs.shape[0]:
                     raise ValueError("Input should be of size (lin_dim,).")
+                xs = reshape(xs, (1, self.lin_dim))
                 n_inputs = 1
             else:
                 n_inputs = xs.shape[0]
diff --git a/tests/distributions/test_abstract_hypercylindrical_distribution.py b/tests/distributions/test_abstract_hypercylindrical_distribution.py
index c375240f4..e4d136499 100644
--- a/tests/distributions/test_abstract_hypercylindrical_distribution.py
+++ b/tests/distributions/test_abstract_hypercylindrical_distribution.py
@@ -65,6 +65,18 @@ def test_get_reasonable_integration_boundaries(self):
             integration_boundaries, array([[0.0, 2.0 * pi], [-4.0, 8.0]])
         )
 
+    def test_integrate_rejects_wrong_boundary_count(self):
+        dist = DummyHypercylindricalDistribution()
+
+        with self.assertRaisesRegex(ValueError, "one .* interval"):
+            dist.integrate_numerically([[0.0, 1.0]])
+
+    def test_integrate_rejects_reversed_boundaries(self):
+        dist = DummyHypercylindricalDistribution()
+
+        with self.assertRaisesRegex(ValueError, "lower bound"):
+            dist.integrate_numerically([[1.0, 0.0], [-1.0, 1.0]])
+
     def test_constructor_accepts_fully_periodic_dimension(self):
         dist = PartiallyWrappedNormalDistribution(array([1.0]), array([[1.0]]), 1)
 
@@ -101,6 +113,36 @@ def test_condition_on_periodic_rejects_wrong_input_dimension(self):
         with self.assertRaisesRegex(ValueError, "bound_dim"):
             hwn.condition_on_periodic(array([1.0]), normalize=False)
 
+    def test_condition_on_linear_accepts_single_multidimensional_periodic_point(self):
+        hwn = PartiallyWrappedNormalDistribution(
+            array([1.0, 2.0, 3.0]),
+            array([[2.0, 0.3, 0.1], [0.3, 1.5, 0.2], [0.1, 0.2, 1.0]]),
+            2,
+        )
+        periodic_point = array([1.5, 2.5])
+
+        dist_cond = hwn.condition_on_linear(array([3.5]), normalize=False)
+
+        npt.assert_allclose(
+            dist_cond.pdf(periodic_point),
+            hwn.pdf(array([periodic_point[0], periodic_point[1], 3.5])),
+        )
+
+    def test_condition_on_periodic_accepts_single_multidimensional_linear_point(self):
+        hwn = PartiallyWrappedNormalDistribution(
+            array([1.0, 2.0, 3.0]),
+            array([[2.0, 0.3, 0.1], [0.3, 1.5, 0.2], [0.1, 0.2, 1.0]]),
+            1,
+        )
+        linear_point = array([2.5, 3.5])
+
+        dist_cond = hwn.condition_on_periodic(array(1.5), normalize=False)
+
+        npt.assert_allclose(
+            dist_cond.pdf(linear_point),
+            hwn.pdf(array([1.5, linear_point[0], linear_point[1]])),
+        )
+
     def test_mode_numerical_rejects_unsupported_backend(self):
         hwn = PartiallyWrappedNormalDistribution(
             array([1.0, 2.0]), array([[2.0, 0.3], [0.3, 1.0]]), 1