We propose and analyse randomized cubature formulae for the numerical integration of functions with respect to a given probability measure μ defined on a domain Γ⊆Rd, in any dimension d. Each cubature formula is exact on a given finite-dimensional subspace Vn⊂L2(Γ,μ) of dimension n, and uses pointwise evaluations of the integrand function ϕ:Γ→R at m>n independent random points. These points are drawn from a suitable auxiliary probability measure that depends on Vn. We show that, up to a logarithmic factor, a linear proportionality between m and n with dimension-independent constant ensures stability of the cubature formula with high probability. We also prove error estimates in probability and in expectation for any n≥1 and m>n, thus covering both preasymptotic and asymptotic regimes. Our analysis shows that the expected cubature error decays as n/m times the L2(Γ,μ)-best approximation error of ϕ in Vn. On the one hand, for fixed n and m→∞ our cubature formula can be seen as a variance reduction technique for a Monte Carlo estimator, and can lead to enormous variance reduction for smooth integrand functions and subspaces Vn with spectral approximation properties. On the other hand, when n,m→∞, our cubature becomes of high order with spectral convergence. As a further contribution, we analyse also another cubature whose expected error decays as 1/m times the L2(Γ,μ)-best approximation error of ϕ in Vn, but with constants that can be larger in the preasymptotic regime. Finally we show that, under a more demanding (at least quadratic) proportionality between m and n, all the weights of the cubature are strictly positive with high probability. As an example of application, we discuss the case where the domain Γ has the structure of Cartesian product, μ is a product measure on Γ and Vn contains algebraic multivariate polynomials.