We propose and analyse randomized cubature formulae for the numerical integration of functions with respect to a given probability measure μ defined on a domain Γ⊆ℝ^d, in any dimension d. Each cubature formula is conceived to be exact on a given finite dimensional subspace V_n⊂L^2(Γ,μ) of dimension n, and uses pointwise evaluations of the integrand function φ:Γ→ℝ at m>n independent random points. These points are distributed according to a suitable auxiliary probability measure that depends on V_n. 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 very high probability. We also prove error estimates in probability and in expectation for any n≥1 and m>n, thus covering both pre-asymptotic and asymptotic regimes. Our analysis shows that the expected cubature error decays as √(n/m) times the L^2(Γ,μ)-best approximation error of φ in V_n. 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 V_n with spectral approximation properties. On the other hand, when we let n,m→∞, our cubature becomes of high order with spectral convergence. Finally we show that, under a more demanding (at least quadratic) proportionality between m and n, the weights of the cubature are positive with very 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 the space V_n contains algebraic multivariate polynomials.