We present cubature formulae for the integration of functions in arbitrary dimension and arbitrary domain. These cubatures are exact on a given finite-dimensional subspace Vn of L^2 of dimension n, they are stable with high probability and are constructed using m pointwise evaluations of the integrand function with m proportional to nlog(n). For these cubatures we provide a convergence analysis showing that the expected cubature error decays as m^{-1/2} times the L^2 best approximation of the integrand function in Vn.