diff --git a/src/probnum/problems/zoo/__init__.py b/src/probnum/problems/zoo/__init__.py
index 951eae1e5..688ebc683 100644
--- a/src/probnum/problems/zoo/__init__.py
+++ b/src/probnum/problems/zoo/__init__.py
@@ -1,7 +1,6 @@
 """Collection of test problems for probabilistic numerical methods.
 
-Subpackage offering implementations of standard example problems or
-convenient interfaces to benchmark problems for probabilistic numerical
-methods. These test problems are meant for rapid experimentation and
-prototyping.
+Subpackage offering implementations of standard example problems or convenient
+interfaces to benchmark problems for probabilistic numerical methods. These test
+problems are meant for rapid experimentation and prototyping.
 """
diff --git a/src/probnum/problems/zoo/quad/__init__.py b/src/probnum/problems/zoo/quad/__init__.py
index 987103016..b5e8e213b 100644
--- a/src/probnum/problems/zoo/quad/__init__.py
+++ b/src/probnum/problems/zoo/quad/__init__.py
@@ -1,7 +1,6 @@
 """Test problems for numerical integration/ quadrature."""
 
 from ._emukit_problems import circulargaussian2d, hennig1d, hennig2d, sombrero2d
-
 from ._quadproblems_gaussian import *
 from ._quadproblems_uniform import *
 
diff --git a/src/probnum/problems/zoo/quad/_quadproblems_gaussian.py b/src/probnum/problems/zoo/quad/_quadproblems_gaussian.py
index 31ab179cf..0d4a764d3 100644
--- a/src/probnum/problems/zoo/quad/_quadproblems_gaussian.py
+++ b/src/probnum/problems/zoo/quad/_quadproblems_gaussian.py
@@ -7,6 +7,7 @@
 from scipy.stats import norm
 
 from probnum.problems import QuadratureProblem
+from probnum.quad.integration_measures import GaussianMeasure
 from probnum.typing import FloatLike
 
 __all__ = [
@@ -65,13 +66,12 @@ def uniform_to_gaussian_quadprob(
        of the 14th Monte Carlo and Quasi-Monte Carlo Methods (MCQMC) conference 2020.
        arXiv:2006.07487.
     """
-
-    dim = len(quadprob.lower_bd)
+    lower_bd, upper_bd = quadprob.measure.domain
 
     # Check that the original quadrature problem was defined on [0,1]^d
-    if np.any(quadprob.lower_bd != 0.0):
+    if np.any(lower_bd != 0.0):
         raise ValueError("quadprob is not an integration problem over [0,1]^d")
-    if np.any(quadprob.upper_bd != 1.0):
+    if np.any(upper_bd != 1.0):
         raise ValueError("quadprob is not an integration problem over [0,1]^d")
 
     # Check that the original quadrature problem has a scalar valued solution
@@ -80,7 +80,7 @@ def uniform_to_gaussian_quadprob(
 
     # Construct transformation of the integrand
     def uniform_to_gaussian_integrand(
-        func: Callable[[np.ndarray], np.ndarray],
+        fun: Callable[[np.ndarray], np.ndarray],
         mean: FloatLike = 0.0,
         var: FloatLike = 1.0,
     ) -> Callable[[np.ndarray], np.ndarray]:
@@ -92,17 +92,14 @@ def uniform_to_gaussian_integrand(
             raise TypeError("The variance should be a positive float.")
 
         def newfunc(x):
-            return func(norm.cdf((x - mean) / var))
+            return fun(norm.cdf((x - mean) / var))
 
         return newfunc
 
+    gaussian_measure = GaussianMeasure(mean=mean, cov=var)
     return QuadratureProblem(
-        integrand=uniform_to_gaussian_integrand(
-            func=quadprob.integrand, mean=mean, var=var
-        ),
-        lower_bd=np.broadcast_to(-np.Inf, dim),
-        upper_bd=np.broadcast_to(np.Inf, dim),
-        output_dim=None,
+        fun=uniform_to_gaussian_integrand(fun=quadprob.fun, mean=mean, var=var),
+        measure=gaussian_measure,
         solution=quadprob.solution,
     )
 
@@ -154,6 +151,18 @@ def sum_polynomials(
     if b is None:
         b = np.broadcast_to(1, (1, dim))
 
+    if len(a.shape) != 2:
+        raise ValueError(
+            f"Invalid shape {a.shape} for parameter `a`. "
+            f"Expected parameters of shape (p+1)xdim"
+        )
+
+    if len(b.shape) != 2:
+        raise ValueError(
+            f"Invalid shape {a.shape} for parameter `a`. "
+            f"Expected parameters of shape (p+1)xdim"
+        )
+
     # Check that the parameters have valid values and shape
     if a.shape[1] != dim:
         raise ValueError(
@@ -187,11 +196,11 @@ 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)
 
+    gaussian_measure = GaussianMeasure(mean=mean, cov=var)
     return QuadratureProblem(
-        integrand=integrand,
-        lower_bd=np.broadcast_to(-np.Inf, dim),
-        upper_bd=np.broadcast_to(np.Inf, dim),
-        output_dim=None,
+        fun=integrand,
+        measure=gaussian_measure,
         solution=solution,
     )
diff --git a/src/probnum/problems/zoo/quad/_quadproblems_uniform.py b/src/probnum/problems/zoo/quad/_quadproblems_uniform.py
index c4f577d3c..250452910 100644
--- a/src/probnum/problems/zoo/quad/_quadproblems_uniform.py
+++ b/src/probnum/problems/zoo/quad/_quadproblems_uniform.py
@@ -21,6 +21,8 @@
     "morokoff_caflisch_2",
 ]
 
+from probnum.quad.integration_measures import LebesgueMeasure
+
 
 def genz_continuous(
     dim: int, a: np.ndarray = None, u: np.ndarray = None
@@ -77,7 +79,7 @@ def genz_continuous(
     if np.any(u < 0.0) or np.any(u > 1):
         raise ValueError("The parameters `u` must lie in the interval [0.0, 1.0].")
 
-    def integrand(x: np.ndarray) -> np.ndarray:
+    def fun(x: np.ndarray) -> np.ndarray:
         nonlocal a, u
         n = x.shape[0]
 
@@ -99,12 +101,12 @@ def integrand(x: np.ndarray) -> np.ndarray:
         return f.reshape((n, 1))
 
     solution = np.prod((2.0 - np.exp(-1.0 * a * u) - np.exp(a * (u - 1))) / a)
-
+    lower_bd = np.broadcast_to(0.0, dim)
+    upper_bd = np.broadcast_to(1.0, dim)
+    uniform_measure = LebesgueMeasure(domain=(lower_bd, upper_bd), input_dim=dim)
     return QuadratureProblem(
-        integrand=integrand,
-        lower_bd=np.broadcast_to(0.0, dim),
-        upper_bd=np.broadcast_to(1.0, dim),
-        output_dim=None,
+        fun=fun,
+        measure=uniform_measure,
         solution=solution,
     )
 
@@ -153,7 +155,7 @@ def genz_cornerpeak(
     if np.any(u < 0.0) or np.any(u > 1):
         raise ValueError("The parameters `u` must lie in the interval [0.0, 1.0].")
 
-    def integrand(x: np.ndarray) -> np.ndarray:
+    def fun(x: np.ndarray) -> np.ndarray:
         nonlocal a, u
         n = x.shape[0]
 
@@ -182,11 +184,12 @@ def integrand(x: np.ndarray) -> np.ndarray:
             ) ** (-1)
     solution = solution / (np.prod(a) * np.math.factorial(dim))
 
+    lower_bd = np.broadcast_to(0.0, dim)
+    upper_bd = np.broadcast_to(1.0, dim)
+    uniform_measure = LebesgueMeasure(domain=(lower_bd, upper_bd), input_dim=dim)
     return QuadratureProblem(
-        integrand=integrand,
-        lower_bd=np.broadcast_to(0.0, dim),
-        upper_bd=np.broadcast_to(1.0, dim),
-        output_dim=None,
+        fun=fun,
+        measure=uniform_measure,
         solution=solution,
     )
 
@@ -240,7 +243,7 @@ def genz_discontinuous(
     if np.any(u < 0.0) or np.any(u > 1):
         raise ValueError("The parameters `u` must lie in the interval [0.0, 1.0].")
 
-    def integrand(x: np.ndarray) -> np.ndarray:
+    def fun(x: np.ndarray) -> np.ndarray:
         nonlocal a, u
         n = x.shape[0]
 
@@ -265,11 +268,12 @@ def integrand(x: np.ndarray) -> np.ndarray:
     if dim > 1:
         solution = np.prod((np.exp(a * np.minimum(u, 1.0)) - 1.0) / a)
 
+    lower_bd = np.broadcast_to(0.0, dim)
+    upper_bd = np.broadcast_to(1.0, dim)
+    uniform_measure = LebesgueMeasure(domain=(lower_bd, upper_bd), input_dim=dim)
     return QuadratureProblem(
-        integrand=integrand,
-        lower_bd=np.broadcast_to(0.0, dim),
-        upper_bd=np.broadcast_to(1.0, dim),
-        output_dim=None,
+        fun=fun,
+        measure=uniform_measure,
         solution=solution,
     )
 
@@ -318,7 +322,7 @@ def genz_gaussian(
     if np.any(u < 0.0) or np.any(u > 1):
         raise ValueError("The parameters `u` must lie in the interval [0.0, 1.0].")
 
-    def integrand(x: np.ndarray) -> np.ndarray:
+    def fun(x: np.ndarray) -> np.ndarray:
         nonlocal a, u
         n = x.shape[0]
 
@@ -343,11 +347,12 @@ def integrand(x: np.ndarray) -> np.ndarray:
         (norm.cdf(np.sqrt(2) * a * (1.0 - u)) - norm.cdf(-np.sqrt(2) * a * u)) / a
     )
 
+    lower_bd = np.broadcast_to(0.0, dim)
+    upper_bd = np.broadcast_to(1.0, dim)
+    uniform_measure = LebesgueMeasure(domain=(lower_bd, upper_bd), input_dim=dim)
     return QuadratureProblem(
-        integrand=integrand,
-        lower_bd=np.broadcast_to(0.0, dim),
-        upper_bd=np.broadcast_to(1.0, dim),
-        output_dim=None,
+        fun=fun,
+        measure=uniform_measure,
         solution=solution,
     )
 
@@ -399,7 +404,7 @@ def genz_oscillatory(
     if np.any(u < 0.0) or np.any(u > 1):
         raise ValueError("The parameters `u` must lie in the interval [0.0, 1.0].")
 
-    def integrand(x: np.ndarray) -> np.ndarray:
+    def fun(x: np.ndarray) -> np.ndarray:
         nonlocal a, u
         n = x.shape[0]
 
@@ -443,11 +448,12 @@ def hfunc(x):
 
     solution = solution / np.prod(a)
 
+    lower_bd = np.broadcast_to(0.0, dim)
+    upper_bd = np.broadcast_to(1.0, dim)
+    uniform_measure = LebesgueMeasure(domain=(lower_bd, upper_bd), input_dim=dim)
     return QuadratureProblem(
-        integrand=integrand,
-        lower_bd=np.broadcast_to(0.0, dim),
-        upper_bd=np.broadcast_to(1.0, dim),
-        output_dim=None,
+        fun=fun,
+        measure=uniform_measure,
         solution=solution,
     )
 
@@ -497,7 +503,7 @@ def genz_productpeak(
     if np.any(u < 0.0) or np.any(u > 1.0):
         raise ValueError("The parameters `u` must lie in the interval [0.0, 1.0].")
 
-    def integrand(x: np.ndarray) -> np.ndarray:
+    def fun(x: np.ndarray) -> np.ndarray:
         nonlocal a, u
         n = x.shape[0]
 
@@ -520,11 +526,12 @@ def integrand(x: np.ndarray) -> np.ndarray:
 
     solution = np.prod(a * (np.arctan(a * (1.0 - u)) - np.arctan(-1.0 * a * u)))
 
+    lower_bd = np.broadcast_to(0.0, dim)
+    upper_bd = np.broadcast_to(1.0, dim)
+    uniform_measure = LebesgueMeasure(domain=(lower_bd, upper_bd), input_dim=dim)
     return QuadratureProblem(
-        integrand=integrand,
-        lower_bd=np.broadcast_to(0.0, dim),
-        upper_bd=np.broadcast_to(1.0, dim),
-        output_dim=None,
+        fun=fun,
+        measure=uniform_measure,
         solution=solution,
     )
 
@@ -549,7 +556,7 @@ def bratley1992(dim: int) -> QuadratureProblem:
     .. [2] https://www.sfu.ca/~ssurjano/bratleyetal92.html
     """
 
-    def integrand(x: np.ndarray) -> np.ndarray:
+    def fun(x: np.ndarray) -> np.ndarray:
         n = x.shape[0]
 
         # Check that the input points have valid values
@@ -571,11 +578,12 @@ def integrand(x: np.ndarray) -> np.ndarray:
 
     solution = -(1.0 / 3) * (1.0 - ((-0.5) ** dim))
 
+    lower_bd = np.broadcast_to(0.0, dim)
+    upper_bd = np.broadcast_to(1.0, dim)
+    uniform_measure = LebesgueMeasure(domain=(lower_bd, upper_bd), input_dim=dim)
     return QuadratureProblem(
-        integrand=integrand,
-        lower_bd=np.broadcast_to(0.0, dim),
-        upper_bd=np.broadcast_to(1.0, dim),
-        output_dim=None,
+        fun=fun,
+        measure=uniform_measure,
         solution=solution,
     )
 
@@ -599,7 +607,7 @@ def roos_arnold(dim: int) -> QuadratureProblem:
     .. [2] https://www.sfu.ca/~ssurjano/roosarn63.html
     """
 
-    def integrand(x: np.ndarray) -> np.ndarray:
+    def fun(x: np.ndarray) -> np.ndarray:
         n = x.shape[0]
 
         # Check that the input points have valid values
@@ -618,11 +626,12 @@ def integrand(x: np.ndarray) -> np.ndarray:
 
     solution = 1.0
 
+    lower_bd = np.broadcast_to(0.0, dim)
+    upper_bd = np.broadcast_to(1.0, dim)
+    uniform_measure = LebesgueMeasure(domain=(lower_bd, upper_bd), input_dim=dim)
     return QuadratureProblem(
-        integrand=integrand,
-        lower_bd=np.broadcast_to(0.0, dim),
-        upper_bd=np.broadcast_to(1.0, dim),
-        output_dim=None,
+        fun=fun,
+        measure=uniform_measure,
         solution=solution,
     )
 
@@ -647,7 +656,7 @@ def gfunction(dim: int) -> QuadratureProblem:
     .. [2] https://www.sfu.ca/~ssurjano/gfunc.html
     """
 
-    def integrand(x: np.ndarray) -> np.ndarray:
+    def fun(x: np.ndarray) -> np.ndarray:
         n = x.shape[0]
 
         # Check that the input points have valid values
@@ -667,11 +676,12 @@ def integrand(x: np.ndarray) -> np.ndarray:
 
     solution = 1.0
 
+    lower_bd = np.broadcast_to(0.0, dim)
+    upper_bd = np.broadcast_to(1.0, dim)
+    uniform_measure = LebesgueMeasure(domain=(lower_bd, upper_bd), input_dim=dim)
     return QuadratureProblem(
-        integrand=integrand,
-        lower_bd=np.broadcast_to(0.0, dim),
-        upper_bd=np.broadcast_to(1.0, dim),
-        output_dim=None,
+        fun=fun,
+        measure=uniform_measure,
         solution=solution,
     )
 
@@ -697,7 +707,7 @@ def morokoff_caflisch_1(dim: int) -> QuadratureProblem:
     .. [3] https://www.sfu.ca/~ssurjano/morcaf95a.html
     """
 
-    def integrand(x: np.ndarray) -> np.ndarray:
+    def fun(x: np.ndarray) -> np.ndarray:
         n = x.shape[0]
 
         # Check that the input points have valid values
@@ -716,11 +726,12 @@ def integrand(x: np.ndarray) -> np.ndarray:
 
     solution = 1.0
 
+    lower_bd = np.broadcast_to(0.0, dim)
+    upper_bd = np.broadcast_to(1.0, dim)
+    uniform_measure = LebesgueMeasure(domain=(lower_bd, upper_bd), input_dim=dim)
     return QuadratureProblem(
-        integrand=integrand,
-        lower_bd=np.broadcast_to(0.0, dim),
-        upper_bd=np.broadcast_to(1.0, dim),
-        output_dim=None,
+        fun=fun,
+        measure=uniform_measure,
         solution=solution,
     )
 
@@ -744,7 +755,7 @@ def morokoff_caflisch_2(dim: int) -> QuadratureProblem:
     .. [2] https://www.sfu.ca/~ssurjano/morcaf95b.html
     """
 
-    def integrand(x: np.ndarray) -> np.ndarray:
+    def fun(x: np.ndarray) -> np.ndarray:
         n = x.shape[0]
 
         # Check that the input points have valid values
@@ -763,10 +774,11 @@ def integrand(x: np.ndarray) -> np.ndarray:
 
     solution = 1.0
 
+    lower_bd = np.broadcast_to(0.0, dim)
+    upper_bd = np.broadcast_to(1.0, dim)
+    uniform_measure = LebesgueMeasure(domain=(lower_bd, upper_bd), input_dim=dim)
     return QuadratureProblem(
-        integrand=integrand,
-        lower_bd=np.broadcast_to(0.0, dim),
-        upper_bd=np.broadcast_to(1.0, dim),
-        output_dim=None,
+        fun=fun,
+        measure=uniform_measure,
         solution=solution,
     )
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 b062b43a7..1d1176634 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,6 +1,9 @@
+from typing import Callable
+
 import numpy as np
 import pytest
 
+from probnum.problems import QuadratureProblem
 from probnum.problems.zoo.quad import (
     bratley1992,
     genz_continuous,
@@ -57,14 +60,16 @@ def test_genz_uniform_param_checks(genz_problem):
         roos_arnold,
     ],
 )
-def test_integrand_eval_checks(quad_problem_constructor):
+def test_integrand_eval_checks(
+    quad_problem_constructor: Callable[..., QuadratureProblem]
+):
     quad_problem = quad_problem_constructor(2)
 
     with pytest.raises(ValueError):
-        quad_problem.integrand(np.zeros((4, 3)))
+        quad_problem.fun(np.zeros((4, 3)))
 
     with pytest.raises(ValueError):
-        quad_problem.integrand(np.full((4, 2), -0.1))
+        quad_problem.fun(np.full((4, 2), -0.1))
 
 
 @pytest.mark.parametrize(
@@ -109,32 +114,32 @@ def test_genz_uniform(a, u, dim):
     # 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.integrand(x_unif)) / n,
+        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.integrand(x_unif)) / n,
+        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.integrand(x_unif)) / n,
+        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.integrand(x_unif)) / n,
+        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.integrand(x_unif)) / n,
+        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.integrand(x_unif)) / n,
+        np.sum(quadprob_genz_productpeak.fun(x_unif)) / n,
         quadprob_genz_productpeak.solution,
         rtol=1e-03,
     )
@@ -167,27 +172,27 @@ def test_integrands_uniform(dim):
     # Test that all Monte Carlo estimators approximate the true value of the integral
     # (integration against uniform)
     np.testing.assert_allclose(
-        np.sum(quadprob_bratley1992.integrand(x_unif)) / n,
+        np.sum(quadprob_bratley1992.fun(x_unif)) / n,
         quadprob_bratley1992.solution,
         rtol=1e-03,
     )
     np.testing.assert_allclose(
-        np.sum(quadprob_roos_arnold.integrand(x_unif)) / n,
+        np.sum(quadprob_roos_arnold.fun(x_unif)) / n,
         quadprob_roos_arnold.solution,
         rtol=1e-03,
     )
     np.testing.assert_allclose(
-        np.sum(quadprob_gfunction.integrand(x_unif)) / n,
+        np.sum(quadprob_gfunction.fun(x_unif)) / n,
         quadprob_gfunction.solution,
         rtol=1e-03,
     )
     np.testing.assert_allclose(
-        np.sum(quadprob_morokoff_caflisch_1.integrand(x_unif)) / n,
+        np.sum(quadprob_morokoff_caflisch_1.fun(x_unif)) / n,
         quadprob_morokoff_caflisch_1.solution,
         rtol=1e-03,
     )
     np.testing.assert_allclose(
-        np.sum(quadprob_morokoff_caflisch_2.integrand(x_unif)) / n,
+        np.sum(quadprob_morokoff_caflisch_2.fun(x_unif)) / n,
         quadprob_morokoff_caflisch_2.solution,
         rtol=1e-03,
     )