Pre-existing methods for multidimensional numerical integration tend to use many more evaluations of the integrand than necessary. That is, pre-existing methods are evaluation inefficient. In some applications of multidimensional numerical integration, such as cosmological model selection, each evaluation is expensive. Consequently, multidimensional numerical integration by way of any pre-existing method can be prohibitively expensive. An evaluationefficient method for multidimensional numerical integration could enable such applications. Presented in this thesis is VoroInt (Voronoi Integrator), which is a novel Monte Carlo method that uses a Voronoi decomposition of the integration domain into regions about sample sites. The approximate integral returned by the method is a Riemann sum of the value at each site weighted by the size (that is, area, volume, or hyper-volume) of the Voronoi region that contains the site. The error estimate returned by the method is the sum of, for each region, the error estimate for the approximate integral over the region. The error estimate for a region is the average deviation of the values at the sites in the regions that are adjacent to the region from the value at the site in the region times the size of the region. A similar formula in which the average deviation is replaced by the absolute maximal deviation (because the latter is more robust than the former) is used by the method to prioritize each region for further sampling. Also presented are comparisons of the number of evaluations used by and the accuracy of the values returned by pre-existing methods and the novel method for various test integrands, which show that the evaluation efficiency and the accuracy of the novel method are better than those of pre-existing methods. The improvements come from the effective sampling and error estimation by the method. Effective sampling enables the method to efficiently find features of the integrand that significantly contribute to the integral. Effective error estimation enables the method to well estimate the accuracy of the approximate integral.