diff --git a/src/probnum/problems/zoo/quad/__init__.py b/src/probnum/problems/zoo/quad/__init__.py
index b5e8e213b..16320cb2c 100644
--- a/src/probnum/problems/zoo/quad/__init__.py
+++ b/src/probnum/problems/zoo/quad/__init__.py
@@ -5,4 +5,22 @@
 from ._quadproblems_uniform import *
 
 # Public classes and functions. Order is reflected in documentation.
-__all__ = ["circulargaussian2d", "hennig1d", "hennig2d", "sombrero2d"]
+__all__ = [
+    "circulargaussian2d",
+    "hennig1d",
+    "hennig2d",
+    "sombrero2d",
+    "uniform_to_gaussian_quadprob",
+    "sum_polynomials",
+    "genz_continuous",
+    "genz_cornerpeak",
+    "genz_discontinuous",
+    "genz_gaussian",
+    "genz_oscillatory",
+    "genz_productpeak",
+    "bratley1992",
+    "roos_arnold",
+    "gfunction",
+    "morokoff_caflisch_1",
+    "morokoff_caflisch_2",
+]
diff --git a/src/probnum/problems/zoo/quad/_quadproblems_gaussian.py b/src/probnum/problems/zoo/quad/_quadproblems_gaussian.py
index 0d4a764d3..aa300a9b3 100644
--- a/src/probnum/problems/zoo/quad/_quadproblems_gaussian.py
+++ b/src/probnum/problems/zoo/quad/_quadproblems_gaussian.py
@@ -1,6 +1,6 @@
 """Test problems for integration against a Gaussian measure."""
 
-from typing import Callable
+from typing import Callable, Union
 
 import numpy as np
 from scipy import special
@@ -16,8 +16,35 @@
 ]
 
 
+# Construct transformation of the integrand
+def uniform_to_gaussian_integrand(
+    fun: Callable[[np.ndarray], np.ndarray],
+    mean: Union[float, np.floating, np.ndarray] = 0.0,
+    std: Union[float, np.floating, np.ndarray] = 1.0,
+) -> Callable[[np.ndarray], np.ndarray]:
+    # mean and var should be either one-dimensional
+    if isinstance(mean, np.ndarray):
+        if len(mean.shape) != 1:
+            raise TypeError(
+                "The mean parameter should be a float or a d-dimensional array."
+            )
+
+    if isinstance(std, np.ndarray):
+        if len(std.shape) != 1:
+            raise TypeError(
+                "The std parameter should be a float or a d-dimensional array."
+            )
+
+    def new_func(x):
+        return fun(norm.cdf(x, loc=mean, scale=std))
+
+    return new_func
+
+
 def uniform_to_gaussian_quadprob(
-    quadprob: QuadratureProblem, mean: FloatLike = 0.0, var: FloatLike = 1.0
+    quadprob: QuadratureProblem,
+    mean: Union[float, np.floating, np.ndarray] = 0.0,
+    std: Union[float, np.floating, np.ndarray] = 1.0,
 ) -> QuadratureProblem:
     r"""Creates a new QuadratureProblem for integration against a Gaussian on
     :math:`\mathbb{R}^d` by using an existing QuadratureProblem whose integrand is
@@ -39,9 +66,12 @@ def uniform_to_gaussian_quadprob(
     quadprob
         A QuadratureProblem instance which includes an integrand defined on [0,1]^d
     mean
-        Mean of the Gaussian distribution.
-    var
-        Diagonal element for the covariance matrix of the Gaussian distribution.
+        Mean of the Gaussian distribution. If `float`, mean is set to the same value
+        across all dimensions. Else, specifies the mean as a d-dimensional array.
+    std
+        Diagonal element for the covariance matrix of the Gaussian distribution. If
+        `float`, the covariance matrix has the same diagonal value for all dimensions.
+        Else, specifies the covariance matrix via a d-dimensional array.
 
     Returns
     -------
@@ -78,27 +108,14 @@ def uniform_to_gaussian_quadprob(
     if np.ndim(quadprob.solution) != 0:
         raise ValueError("The solution of quadprob is not a scalar.")
 
-    # Construct transformation of the integrand
-    def uniform_to_gaussian_integrand(
-        fun: Callable[[np.ndarray], np.ndarray],
-        mean: FloatLike = 0.0,
-        var: FloatLike = 1.0,
-    ) -> Callable[[np.ndarray], np.ndarray]:
-        # mean and var should be either one-dimensional, or an array of dimension d
-        if isinstance(mean, float) is False:
-            raise TypeError("The mean parameter should be a float.")
-
-        if isinstance(var, float) is False or var <= 0.0:
-            raise TypeError("The variance should be a positive float.")
-
-        def newfunc(x):
-            return fun(norm.cdf((x - mean) / var))
+    dim = lower_bd.shape[0]
 
-        return newfunc
-
-    gaussian_measure = GaussianMeasure(mean=mean, cov=var)
+    cov = std**2
+    if isinstance(std, np.ndarray):
+        cov = np.eye(dim) * cov
+    gaussian_measure = GaussianMeasure(mean=mean, cov=cov, input_dim=dim)
     return QuadratureProblem(
-        fun=uniform_to_gaussian_integrand(fun=quadprob.fun, mean=mean, var=var),
+        fun=uniform_to_gaussian_integrand(fun=quadprob.fun, mean=mean, std=std),
         measure=gaussian_measure,
         solution=quadprob.solution,
     )
@@ -159,7 +176,7 @@ def sum_polynomials(
 
     if len(b.shape) != 2:
         raise ValueError(
-            f"Invalid shape {a.shape} for parameter `a`. "
+            f"Invalid shape {b.shape} for parameter `b`. "
             f"Expected parameters of shape (p+1)xdim"
         )
 
@@ -196,9 +213,12 @@ def integrand(x: np.ndarray) -> np.ndarray:
     delta = (np.remainder(b, 2) - 1) ** 2
     doublefact = special.factorial2(b - 1)
     solution = np.sum(np.prod(a * delta * (var**b) * doublefact, axis=1))
-    mean = np.zeros_like(dim)
+    if isinstance(var, float):
+        mean = 0.0
+    else:
+        mean = np.zeros(dim)
 
-    gaussian_measure = GaussianMeasure(mean=mean, cov=var)
+    gaussian_measure = GaussianMeasure(mean=mean, cov=var, input_dim=dim)
     return QuadratureProblem(
         fun=integrand,
         measure=gaussian_measure,
diff --git a/src/probnum/problems/zoo/quad/_quadproblems_uniform.py b/src/probnum/problems/zoo/quad/_quadproblems_uniform.py
index 250452910..7883e49d5 100644
--- a/src/probnum/problems/zoo/quad/_quadproblems_uniform.py
+++ b/src/probnum/problems/zoo/quad/_quadproblems_uniform.py
@@ -264,7 +264,7 @@ def fun(x: np.ndarray) -> np.ndarray:
         return f.reshape(n, 1)
 
     if dim == 1:
-        solution = (np.exp(a * u) - 1.0) / a
+        solution = ((np.exp(a * u) - 1.0) / a).item()
     if dim > 1:
         solution = np.prod((np.exp(a * np.minimum(u, 1.0)) - 1.0) / a)
 
diff --git a/tests/test_problems/test_zoo/test_quad/test_quadproblems_gaussian.py b/tests/test_problems/test_zoo/test_quad/test_quadproblems_gaussian.py
new file mode 100644
index 000000000..df932036b
--- /dev/null
+++ b/tests/test_problems/test_zoo/test_quad/test_quadproblems_gaussian.py
@@ -0,0 +1,239 @@
+from itertools import product
+from typing import Callable, Optional, Union
+
+import numpy as np
+import pytest
+from pytest_cases import parametrize_with_cases
+
+from probnum.problems import QuadratureProblem
+from probnum.problems.zoo.quad import (
+    bratley1992,
+    genz_continuous,
+    genz_cornerpeak,
+    genz_discontinuous,
+    genz_gaussian,
+    genz_oscillatory,
+    genz_productpeak,
+    gfunction,
+    morokoff_caflisch_1,
+    morokoff_caflisch_2,
+    roos_arnold,
+    sum_polynomials,
+    uniform_to_gaussian_quadprob,
+)
+from probnum.quad.integration_measures import LebesgueMeasure
+from probnum.quad.typing import DomainLike
+from tests.test_problems.test_zoo.test_quad.test_quadproblems_gaussian_cases import (
+    GenzStandardNormalCases,
+    GenzVariedNormalCases,
+    OtherIntegrandsGaussianCases,
+)
+
+dim_values = [1, 2]
+quad_prob_constructor_values = [
+    bratley1992,
+    genz_continuous,
+    genz_cornerpeak,
+    genz_discontinuous,
+    genz_gaussian,
+    genz_oscillatory,
+    genz_productpeak,
+    gfunction,
+    morokoff_caflisch_1,
+    morokoff_caflisch_2,
+    roos_arnold,
+]
+param_combinations = list(product(dim_values, quad_prob_constructor_values))
+
+
+@pytest.mark.parametrize("quad_prob_constructor", quad_prob_constructor_values)
+@pytest.mark.parametrize("dim", dim_values)
+def test_wrapping_all_test_functions_works(
+    dim: int,
+    quad_prob_constructor: Callable[[int], QuadratureProblem],
+):
+    """Integration test that wrapping all problems works."""
+    quadprob = quad_prob_constructor(dim)
+    gaussian_quadprob = uniform_to_gaussian_quadprob(quadprob)
+    assert isinstance(gaussian_quadprob, QuadratureProblem)
+
+
+@pytest.mark.parametrize(
+    "domain, dim",
+    [
+        ((0.5, 1.0), 1),
+        ((0.5, 1.0), 2),
+        ((0.0, 1.5), 1),
+        ((0.0, 1.5), 2),
+    ],
+)
+def test_wrong_measure_bounds_scalar_raises(domain: DomainLike, dim: int):
+    uniform_measure = LebesgueMeasure(domain=domain, input_dim=dim)
+    quadprob = genz_continuous(dim)
+    quadprob.measure = uniform_measure
+    with pytest.raises(ValueError) as exc_info:
+        uniform_to_gaussian_quadprob(quadprob)
+
+    assert "[0,1]" in str(exc_info.value)
+
+
+@pytest.mark.parametrize(
+    "domain",
+    [
+        (np.array([0.0, 0.5]), np.array([1.0, 1.0])),
+        (np.array([0.0, 0.0]), np.array([1.0, 1.5])),
+    ],
+)
+def test_wrong_measure_bounds_array_raises(domain: DomainLike):
+    lower_bd, _ = domain
+    dim = len(lower_bd)
+    uniform_measure = LebesgueMeasure(domain=domain, input_dim=dim)
+    quadprob = genz_continuous(dim)
+    quadprob.measure = uniform_measure
+    with pytest.raises(ValueError) as exc_info:
+        uniform_to_gaussian_quadprob(quadprob)
+
+    assert "[0,1]" in str(exc_info.value)
+
+
+def test_wrong_data_shapes_for_mean_raises():
+    with pytest.raises(TypeError) as exc_info:
+        uniform_to_gaussian_quadprob(
+            genz_continuous(2), mean=np.array([[0.0, 0.5]]), std=np.array([0.0, 0.5])
+        )
+
+    assert "mean parameter" in str(exc_info.value)
+
+
+def test_wrong_data_shapes_for_std_raises():
+    with pytest.raises(TypeError) as exc_info:
+        uniform_to_gaussian_quadprob(
+            genz_continuous(2), mean=np.array([0.0, 0.5]), std=np.array([[0.0, 0.5]])
+        )
+
+    assert "std parameter" in str(exc_info.value)
+
+
+def test_wrong_shaped_a_in_sum_polynomials_raises():
+    a = np.array([1.0, 1.0])
+    dim = 2
+    with pytest.raises(ValueError) as exc_info:
+        sum_polynomials(dim=dim, a=a)
+
+    assert "Invalid shape" in str(exc_info.value)
+    assert "parameter `a`" in str(exc_info.value)
+
+
+def test_wrong_shaped_b_in_sum_polynomials_raises():
+    b = np.array([1.0, 1.0])
+    dim = 2
+    with pytest.raises(ValueError) as exc_info:
+        sum_polynomials(dim=dim, b=b)
+
+    assert "Invalid shape" in str(exc_info.value)
+    assert "parameter `b`" in str(exc_info.value)
+
+
+def test_mismatch_dim_a_in_sum_polynomials_raises():
+    a = np.array([[1.0, 1.0]])
+    dim = 1
+    with pytest.raises(ValueError) as exc_info:
+        sum_polynomials(dim=dim, a=a)
+
+    assert "parameter `a`" in str(exc_info.value)
+    assert f"Expected {dim} columns" in str(exc_info.value)
+
+
+def test_mismatch_dim_b_in_sum_polynomials_raises():
+    b = np.array([[1.0, 1.0]])
+    dim = 1
+    with pytest.raises(ValueError) as exc_info:
+        sum_polynomials(dim=dim, b=b)
+
+    assert "parameter `b`" in str(exc_info.value)
+    assert f"Expected {dim} columns" in str(exc_info.value)
+
+
+def test_negative_values_for_b_in_sum_polynomials_raises():
+    b = np.array([[-0.5, 1.0]])
+    dim = 2
+    with pytest.raises(ValueError) as exc_info:
+        sum_polynomials(dim=dim, b=b)
+
+    assert "negative" in str(exc_info.value)
+    assert "parameters `b`" in str(exc_info.value)
+
+
+@parametrize_with_cases(
+    "quadprob, rtol",
+    cases=(
+        GenzStandardNormalCases,
+        GenzVariedNormalCases,
+        OtherIntegrandsGaussianCases,
+    ),
+)
+def test_genz_gaussian(quadprob: QuadratureProblem, rtol: float):
+    """Compare a Monte Carlo estimator against Gaussian measure with a very large number
+    of samples to the true value of the integral.
+
+    The former should be an approximation of the latter.
+    """
+    # Number of Monte Carlo samples for the test
+    n = 10000000
+
+    # generate some normally distributed points from N(0, 1)
+    rng = np.random.default_rng(0)
+    x_gaussian = quadprob.measure.sample(n, rng=rng)  # Monte Carlo samples x_1,...,x_n
+
+    # Test that all Monte Carlo estimators approximate the true value of the integral
+    # (integration against Gaussian)
+    np.testing.assert_allclose(
+        np.sum(quadprob.fun(x_gaussian)) / n,
+        quadprob.solution,
+        rtol=rtol,
+    )
+
+
+@pytest.mark.parametrize(
+    "dim, a, b, var",
+    [
+        (1, None, None, 1.0),
+        (2, None, None, 1.0),
+        (1, np.array([[0.5]]), np.array([[1]]), 0.5),
+        (1, np.array([[0.5]]), np.array([[1]]), np.array([[0.5]])),
+        (
+            2,
+            np.array([[0.5, 1.0], [2.0, 2.0]]),
+            np.array([[1, 2], [2, 3]]),
+            np.array([[1.0, 0], [0, 0.5]]),
+        ),
+    ],
+)
+def test_sum_polynomials(
+    dim: int,
+    a: Optional[np.ndarray],
+    b: Optional[np.ndarray],
+    var: Union[float, np.ndarray],
+):
+    """Compare a Monte Carlo estimator against Gaussian measure with a very large number
+    of samples to the true value of the integral.
+
+    The former should be an approximation of the latter.
+    """
+
+    # Number of Monte Carlo samples for the test
+    n = 10000000
+
+    quadprob = sum_polynomials(dim, a, b, var)
+
+    # generate some normally distributed points from N(0, 1)
+    rng = np.random.default_rng(0)
+    x_gaussian = quadprob.measure.sample(n, rng=rng)  # Monte Carlo samples x_1,...,x_n
+
+    # Test that all Monte Carlo estimators approximate the true value of the integral
+    # (integration against uniform)
+    np.testing.assert_allclose(
+        np.sum(quadprob.fun(x_gaussian)) / n,
+        quadprob.solution,
+        atol=1e-02,
+    )
diff --git a/tests/test_problems/test_zoo/test_quad/test_quadproblems_gaussian_cases.py b/tests/test_problems/test_zoo/test_quad/test_quadproblems_gaussian_cases.py
new file mode 100644
index 000000000..4061684d4
--- /dev/null
+++ b/tests/test_problems/test_zoo/test_quad/test_quadproblems_gaussian_cases.py
@@ -0,0 +1,140 @@
+import numpy as np
+import pytest
+
+from probnum.problems.zoo.quad import (
+    bratley1992,
+    genz_continuous,
+    genz_cornerpeak,
+    genz_discontinuous,
+    genz_gaussian,
+    genz_oscillatory,
+    genz_productpeak,
+    gfunction,
+    morokoff_caflisch_1,
+    roos_arnold,
+    uniform_to_gaussian_quadprob,
+)
+
+genz_params = [
+    (None, None, 3),
+    (5, 0.5, 1),
+    (5.0, 1.0, 1),
+    (2.0, 0.8, 1),
+    (5, 0.5, 4),
+    (5.0, 1.0, 4),
+    (2.0, 0.8, 4),
+]
+
+
+class GenzStandardNormalCases:
+    @pytest.mark.parametrize("a, u, dim", genz_params)
+    def case_genz_continuous(self, dim, a, u):
+        a_vec = np.repeat(a, dim) if a is not None else None
+        u_vec = np.repeat(u, dim) if u is not None else None
+        quadprob = uniform_to_gaussian_quadprob(
+            genz_continuous(dim=dim, a=a_vec, u=u_vec)
+        )
+        rtol = 1e-2
+        return quadprob, rtol
+
+    @pytest.mark.parametrize("a, u, dim", genz_params)
+    def case_genz_cornerpeak(self, dim, a, u):
+        a_vec = np.repeat(a, dim) if a is not None else None
+        u_vec = np.repeat(u, dim) if u is not None else None
+        quadprob = uniform_to_gaussian_quadprob(
+            genz_cornerpeak(dim=dim, a=a_vec, u=u_vec)
+        )
+        rtol = 1e-2
+        return quadprob, rtol
+
+    @pytest.mark.parametrize("a, u, dim", genz_params)
+    def case_genz_discontinuous(self, dim, a, u):
+        a_vec = np.repeat(a, dim) if a is not None else None
+        u_vec = np.repeat(u, dim) if u is not None else None
+        quadprob = uniform_to_gaussian_quadprob(
+            genz_discontinuous(dim=dim, a=a_vec, u=u_vec)
+        )
+        rtol = 1e-2
+        return quadprob, rtol
+
+    @pytest.mark.parametrize("a, u, dim", genz_params)
+    def case_genz_gaussian(self, dim, a, u):
+        a_vec = np.repeat(a, dim) if a is not None else None
+        u_vec = np.repeat(u, dim) if u is not None else None
+        quadprob = uniform_to_gaussian_quadprob(
+            genz_gaussian(dim=dim, a=a_vec, u=u_vec)
+        )
+        rtol = 1e-2
+        return quadprob, rtol
+
+    @pytest.mark.parametrize("a, u, dim", genz_params)
+    def case_genz_oscillatory(self, dim, a, u):
+        a_vec = np.repeat(a, dim) if a is not None else None
+        u_vec = np.repeat(u, dim) if u is not None else None
+        quadprob = uniform_to_gaussian_quadprob(
+            genz_oscillatory(dim=dim, a=a_vec, u=u_vec)
+        )
+        rtol = 9e-2
+        return quadprob, rtol
+
+    @pytest.mark.parametrize("a, u, dim", genz_params)
+    def case_genz_productpeak(self, dim, a, u):
+        a_vec = np.repeat(a, dim) if a is not None else None
+        u_vec = np.repeat(u, dim) if u is not None else None
+        quadprob = uniform_to_gaussian_quadprob(
+            genz_productpeak(dim=dim, a=a_vec, u=u_vec)
+        )
+        rtol = 1e-2
+        return quadprob, rtol
+
+
+normal_params = [
+    (3.0, 2.0, 2),
+    (0.5, 2.0, 1),
+    (np.array([0.0, 1.0]), np.array([0.5, 2.0]), 2),
+]
+
+
+class GenzVariedNormalCases:
+    @pytest.mark.parametrize("mean, std, dim", normal_params)
+    def case_genz_continuous(self, mean, std, dim):
+        quadprob = uniform_to_gaussian_quadprob(
+            genz_continuous(dim=dim), mean=mean, std=std
+        )
+        rtol = 1e-2
+        return quadprob, rtol
+
+
+dim_params = [1, 2]
+
+
+class OtherIntegrandsGaussianCases:
+    @pytest.mark.parametrize("dim", dim_params)
+    def case_bratley1992(self, dim):
+        quadprob = uniform_to_gaussian_quadprob(bratley1992(dim=dim))
+        rtol = 1e-3
+        return quadprob, rtol
+
+    @pytest.mark.parametrize("dim", dim_params)
+    def case_roos_arnold(self, dim):
+        quadprob = uniform_to_gaussian_quadprob(roos_arnold(dim=dim))
+        rtol = 1e-3
+        return quadprob, rtol
+
+    @pytest.mark.parametrize("dim", dim_params)
+    def case_gfunction(self, dim):
+        quadprob = uniform_to_gaussian_quadprob(gfunction(dim=dim))
+        rtol = 1e-3
+        return quadprob, rtol
+
+    @pytest.mark.parametrize("dim", dim_params)
+    def case_morokoff_caflisch_1(self, dim):
+        quadprob = uniform_to_gaussian_quadprob(morokoff_caflisch_1(dim=dim))
+        rtol = 1e-3
+        return quadprob, rtol
+
+    @pytest.mark.parametrize("dim", dim_params)
+    def case_morokoff_caflisch_2(self, dim):
+        quadprob = uniform_to_gaussian_quadprob(morokoff_caflisch_1(dim=dim))
+        rtol = 1e-3
+        return quadprob, rtol
diff --git a/tests/test_problems/test_zoo/test_quad/test_quadproblems_uniform.py b/tests/test_problems/test_zoo/test_quad/test_quadproblems_uniform.py
index 1d1176634..424959796 100644
--- a/tests/test_problems/test_zoo/test_quad/test_quadproblems_uniform.py
+++ b/tests/test_problems/test_zoo/test_quad/test_quadproblems_uniform.py
@@ -1,7 +1,9 @@
+from itertools import product
 from typing import Callable
 
 import numpy as np
 import pytest
+from pytest_cases import parametrize_with_cases
 
 from probnum.problems import QuadratureProblem
 from probnum.problems.zoo.quad import (
@@ -17,6 +19,10 @@
     morokoff_caflisch_2,
     roos_arnold,
 )
+from tests.test_problems.test_zoo.test_quad.test_quadproblems_uniform_cases import (
+    GenzUniformCases,
+    OtherIntegrandsUniformCases,
+)
 
 
 @pytest.mark.parametrize(
@@ -72,83 +78,39 @@ def test_integrand_eval_checks(
         quad_problem.fun(np.full((4, 2), -0.1))
 
 
-@pytest.mark.parametrize(
-    "a, u, dim",
-    [
-        (None, None, 3),
-        (5, 0.5, 1),
-        (5.0, 1.0, 1),
-        (2.0, 0.8, 1),
-        (5, 0.5, 4),
-        (5.0, 1.0, 4),
-        (2.0, 0.8, 4),
-    ],
-)
-def test_genz_uniform(a, u, dim):
-    """Compare a Monte Carlo estimator with a very large number of samples to the true
-    value of the integral.
-
-    The former should be an approximation of the latter.
-    """
-
-    # Number of Monte Carlo samples for the test
-    n = 10000000
-
-    # Define parameters a and u
-    a_vec = np.repeat(a, dim) if a is not None else None
-    u_vec = np.repeat(u, dim) if u is not None else None
-
-    quadprob_genz_continuous = genz_continuous(dim=dim, a=a_vec, u=u_vec)
-    quadprob_genz_cornerpeak = genz_cornerpeak(dim=dim, a=a_vec, u=u_vec)
-    quadprob_genz_discontinuous = genz_discontinuous(dim=dim, a=a_vec, u=u_vec)
-    quadprob_genz_gaussian = genz_gaussian(dim=dim, a=a_vec, u=u_vec)
-    quadprob_genz_oscillatory = genz_oscillatory(dim=dim, a=a_vec, u=u_vec)
-    quadprob_genz_productpeak = genz_productpeak(dim=dim, a=a_vec, u=u_vec)
+dim_values = [1, 2]
+quad_prob_constructor_values = [
+    bratley1992,
+    genz_continuous,
+    genz_cornerpeak,
+    genz_discontinuous,
+    genz_gaussian,
+    genz_oscillatory,
+    genz_productpeak,
+    gfunction,
+    morokoff_caflisch_1,
+    morokoff_caflisch_2,
+    roos_arnold,
+]
+param_combinations = list(product(dim_values, quad_prob_constructor_values))
 
-    # generate some uniformly distributed points on [0,1]^d
-    np.random.seed(0)
-    x_unif = np.random.uniform(
-        low=0.0, high=1.0, size=(n, dim)
-    )  # Monte Carlo samples x_1,...,x_n
 
-    # Test that all Monte Carlo estimators approximate the true value of the integral
-    # (integration against uniform)
-    np.testing.assert_allclose(
-        np.sum(quadprob_genz_continuous.fun(x_unif)) / n,
-        quadprob_genz_continuous.solution,
-        rtol=1e-03,
-    )
-    np.testing.assert_allclose(
-        np.sum(quadprob_genz_cornerpeak.fun(x_unif)) / n,
-        quadprob_genz_cornerpeak.solution,
-        rtol=3e-03,
-    )
-    np.testing.assert_allclose(
-        np.sum(quadprob_genz_discontinuous.fun(x_unif)) / n,
-        quadprob_genz_discontinuous.solution,
-        rtol=1e-03,
-    )
-    np.testing.assert_allclose(
-        np.sum(quadprob_genz_gaussian.fun(x_unif)) / n,
-        quadprob_genz_gaussian.solution,
-        rtol=2e-03,
-    )
-    np.testing.assert_allclose(
-        np.sum(quadprob_genz_oscillatory.fun(x_unif)) / n,
-        quadprob_genz_oscillatory.solution,
-        rtol=9e-02,
-    )
-    np.testing.assert_allclose(
-        np.sum(quadprob_genz_productpeak.fun(x_unif)) / n,
-        quadprob_genz_productpeak.solution,
-        rtol=1e-03,
-    )
+@pytest.mark.parametrize("dim, quad_prob_constructor", param_combinations)
+def test_integrand_solution_float(
+    dim: int, quad_prob_constructor: Callable[..., QuadratureProblem]
+):
+    quadprob = quad_prob_constructor(dim)
+    if np.ndim(quadprob.solution) != 0:
+        raise ValueError(f"The solution of {quadprob} is not a scalar.")
 
 
-@pytest.mark.parametrize("dim", [(1), (4), (10)])
-def test_integrands_uniform(dim):
-    """Compare a Monte Carlo estimator with a very large number of samples to the true
-    value of the integral.
+@parametrize_with_cases(
+    "quadprob, rtol",
+    cases=(GenzUniformCases, OtherIntegrandsUniformCases),
+)
+def test_quadprob_uniform_with_mc(quadprob, rtol):
+    """Compare a Monte Carlo estimator against uniform measure with a very large number
+    of samples to the true value of the integral.
 
     The former should be an approximation of the latter.
     """
@@ -156,43 +118,14 @@ def test_integrands_uniform(dim):
     # Number of Monte Carlo samples for the test
     n = 10000000
 
-    # generate some uniformly distributed points on [0,1]^d
-    np.random.seed(0)
-    x_unif = np.random.uniform(
-        low=0.0, high=1.0, size=(n, dim)
-    )  # Monte Carlo samples x_1,...,x_n
-
-    # Set the integrands to test
-    quadprob_bratley1992 = bratley1992(dim=dim)
-    quadprob_roos_arnold = roos_arnold(dim=dim)
-    quadprob_gfunction = gfunction(dim=dim)
-    quadprob_morokoff_caflisch_1 = morokoff_caflisch_1(dim=dim)
-    quadprob_morokoff_caflisch_2 = morokoff_caflisch_2(dim=dim)
+    # generate some uniformly distributed points from [0, 1]
+    rng = np.random.default_rng(0)
+    x_uniform = quadprob.measure.sample(n, rng=rng)  # Monte Carlo samples x_1,...,x_n
 
     # Test that all Monte Carlo estimators approximate the true value of the integral
-    # (integration against uniform)
-    np.testing.assert_allclose(
-        np.sum(quadprob_bratley1992.fun(x_unif)) / n,
-        quadprob_bratley1992.solution,
-        rtol=1e-03,
-    )
-    np.testing.assert_allclose(
-        np.sum(quadprob_roos_arnold.fun(x_unif)) / n,
-        quadprob_roos_arnold.solution,
-        rtol=1e-03,
-    )
-    np.testing.assert_allclose(
-        np.sum(quadprob_gfunction.fun(x_unif)) / n,
-        quadprob_gfunction.solution,
-        rtol=1e-03,
-    )
-    np.testing.assert_allclose(
-        np.sum(quadprob_morokoff_caflisch_1.fun(x_unif)) / n,
-        quadprob_morokoff_caflisch_1.solution,
-        rtol=1e-03,
-    )
+    # (integration against [0, 1])
     np.testing.assert_allclose(
-        np.sum(quadprob_morokoff_caflisch_2.fun(x_unif)) / n,
-        quadprob_morokoff_caflisch_2.solution,
-        rtol=1e-03,
+        np.sum(quadprob.fun(x_uniform)) / n,
+        quadprob.solution,
+        rtol=rtol,
     )
diff --git a/tests/test_problems/test_zoo/test_quad/test_quadproblems_uniform_cases.py b/tests/test_problems/test_zoo/test_quad/test_quadproblems_uniform_cases.py
new file mode 100644
index 000000000..66f067e10
--- /dev/null
+++ b/tests/test_problems/test_zoo/test_quad/test_quadproblems_uniform_cases.py
@@ -0,0 +1,110 @@
+import numpy as np
+import pytest
+
+from probnum.problems.zoo.quad import (
+    bratley1992,
+    genz_continuous,
+    genz_cornerpeak,
+    genz_discontinuous,
+    genz_gaussian,
+    genz_oscillatory,
+    genz_productpeak,
+    gfunction,
+    morokoff_caflisch_1,
+    roos_arnold, morokoff_caflisch_2,
+)
+
+genz_params = [
+    (None, None, 3),
+    (5, 0.5, 1),
+    (5.0, 1.0, 1),
+    (2.0, 0.8, 1),
+    (5, 0.5, 4),
+    (5.0, 1.0, 4),
+    (2.0, 0.8, 4),
+]
+
+
+class GenzUniformCases:
+    @pytest.mark.parametrize("a, u, dim", genz_params)
+    def case_genz_continuous(self, dim, a, u):
+        a_vec = np.repeat(a, dim) if a is not None else None
+        u_vec = np.repeat(u, dim) if u is not None else None
+        quadprob = genz_continuous(dim=dim, a=a_vec, u=u_vec)
+        rtol = 1e-3
+        return quadprob, rtol
+
+    @pytest.mark.parametrize("a, u, dim", genz_params)
+    def case_genz_cornerpeak(self, dim, a, u):
+        a_vec = np.repeat(a, dim) if a is not None else None
+        u_vec = np.repeat(u, dim) if u is not None else None
+        quadprob = genz_cornerpeak(dim=dim, a=a_vec, u=u_vec)
+        rtol = 1e-2
+        return quadprob, rtol
+
+    @pytest.mark.parametrize("a, u, dim", genz_params)
+    def case_genz_discontinuous(self, dim, a, u):
+        a_vec = np.repeat(a, dim) if a is not None else None
+        u_vec = np.repeat(u, dim) if u is not None else None
+        quadprob = genz_discontinuous(dim=dim, a=a_vec, u=u_vec)
+        rtol = 2e-3
+        return quadprob, rtol
+
+    @pytest.mark.parametrize("a, u, dim", genz_params)
+    def case_genz_gaussian(self, dim, a, u):
+        a_vec = np.repeat(a, dim) if a is not None else None
+        u_vec = np.repeat(u, dim) if u is not None else None
+        quadprob = genz_gaussian(dim=dim, a=a_vec, u=u_vec)
+        rtol = 1e-2
+        return quadprob, rtol
+
+    @pytest.mark.parametrize("a, u, dim", genz_params)
+    def case_genz_oscillatory(self, dim, a, u):
+        a_vec = np.repeat(a, dim) if a is not None else None
+        u_vec = np.repeat(u, dim) if u is not None else None
+        quadprob = genz_oscillatory(dim=dim, a=a_vec, u=u_vec)
+        rtol = 9e-2
+        return quadprob, rtol
+
+    @pytest.mark.parametrize("a, u, dim", genz_params)
+    def case_genz_productpeak(self, dim, a, u):
+        a_vec = np.repeat(a, dim) if a is not None else None
+        u_vec = np.repeat(u, dim) if u is not None else None
+        quadprob = genz_productpeak(dim=dim, a=a_vec, u=u_vec)
+        rtol = 1e-3
+        return quadprob, rtol
+
+
+dim_params = [1, 2, 10]
+
+
+class OtherIntegrandsUniformCases:
+    @pytest.mark.parametrize("dim", dim_params)
+    def case_bratley1992(self, dim):
+        quadprob = bratley1992(dim=dim)
+        rtol = 1e-3
+        return quadprob, rtol
+
+    @pytest.mark.parametrize("dim", dim_params)
+    def case_roos_arnold(self, dim):
+        quadprob = roos_arnold(dim=dim)
+        rtol = 1e-3
+        return quadprob, rtol
+
+    @pytest.mark.parametrize("dim", dim_params)
+    def case_gfunction(self, dim):
+        quadprob = gfunction(dim=dim)
+        rtol = 1e-3
+        return quadprob, rtol
+
+    @pytest.mark.parametrize("dim", dim_params)
+    def case_morokoff_caflisch_1(self, dim):
+        quadprob = morokoff_caflisch_1(dim=dim)
+        rtol = 1e-3
+        return quadprob, rtol
+
+    @pytest.mark.parametrize("dim", dim_params)
+    def case_morokoff_caflisch_2(self, dim):
+        quadprob = morokoff_caflisch_2(dim=dim)
+        rtol = 1e-3
+        return quadprob, rtol