Classical algorithms for numerical integration (quadrature/cubature) proceed by approximating the integrand with a simple function (e.g. a polynomial), and integrate the approximant exactly. In high-dimensional integration, such methods quickly become infeasible due to the curse of dimensionality. A common alternative is the Monte Carlo method (MC), which simply takes the average of random samples, improving the estimate as more and more samples are taken. The main issue with MC is its slow "sqrt(variance/#samples)" convergence, and various techniques have been proposed to reduce the variance. In this work we reveal a numerical analyst's interpretation of MC: it approximates the integrand with a constant function, and integrates the constant exactly. This observation leads naturally to MC-like methods that combines MC with (high-dimensional) function approximation theory, including polynomial approximation, sparse grids and low-rank approximation. The resulting method can be regarded as another variance reduction technique for Monte Carlo.