  We study the complexity of approximating integrals of smooth functions at
absolute precision $\varepsilon > 0$ with confidence level $1 - \delta \in
(0,1)$. The optimal error rate for multivariate functions from classical
isotropic Sobolev spaces $W_p^r(G)$ with sufficient smoothness on bounded
Lipschitz domains $G \subset \mathbb{R}^d$ is determined. It turns out that the
integrability index $p$ has an effect on the influence of the uncertainty
$\delta$ in the complexity. In the limiting case $p = 1$ we see that
deterministic methods cannot be improved by randomization. In general, higher
smoothness reduces the additional effort for diminishing the uncertainty.
Finally, we add a discussion about this problem for function spaces with mixed
smoothness.
