We study the numerical integration of functions from isotropic Sobolev spaces W p s ( [ 0 , 1 ] d ) using finitely many function evaluations within randomized algorithms, aiming for the smallest possible probabilistic error guarantee ε > 0 at confidence level 1 − δ ∈ ( 0 , 1 ) . For spaces consisting of continuous functions, non-linear Monte Carlo methods with optimal confidence properties have already been known, in few cases even linear methods that succeed in that respect. In this paper we promote a method called stratified control variates (SCV) and by it show that already linear methods achieve optimal probabilistic error rates in the high smoothness regime without the need to adjust algorithmic parameters to the uncertainty δ . We also analyse a version of SCV in the low smoothness regime where W p s ( [ 0 , 1 ] d ) may contain functions with singularities. Here, we observe a polynomial dependence of the error on δ − 1 in contrast to the logarithmic dependence in the high smoothness regime.