A Note on Monte Carlo Integration in High Dimensions

Abstract

Monte Carlo integration is a commonly used technique to compute intractable integrals and is typically thought to perform poorly for very high-dimensional integrals. To show that this is not always the case, we examine Monte Carlo integration using techniques from the high-dimensional statistics literature by allowing the dimension of the integral to increase. In doing so, we derive nonasymptotic bounds for the relative and absolute error of the approximation for some general classes of functions through concentration inequalities. We provide concrete examples in which the magnitude of the number of points sampled needed to guarantee a consistent estimate varies between polynomial to exponential, and show that in theory arbitrarily fast or slow rates are possible. This demonstrates that the behavior of Monte Carlo integration in high dimensions is not uniform. Through our methods we also obtain nonasymptotic confidence intervals which are valid regardless of the number of points sampled.

KEYWORDS:

• Concentration inequalities
• High-dimensional statistics
• Numerical integration

1 Introduction

Monte Carlo integration is a simple and widely used numerical integration method which leverages randomness and computational power to approximate intractable integrals. Chapters 1 and 2 of Owen (2013) provide some historical details as well as an introduction to the subject. Roughly speaking, we can think of this method as generating a random grid of points on which we evaluate the function being integrated, and average these values to form an estimate of the integral. As a result, Monte Carlo integration is generally thought to perform poorly for very high-dimensional functions. The intuition behind this is clear: the curse of dimensionality makes covering the region of integration exponentially more difficult as the dimension of the integral increases. However, this poor scaling of the accuracy of the method is not immediate from the usual statement about the error rates of Monte Carlo methods. While it is true that the approximation error is $O(n^{-1/2})$ (Owen 2013, chap. 2.1) on the absolute scale regardless of the dimension of the integrand, where n is the sample size, the constant involved in this rate can depend on the dimension of the integral. Furthermore, for many statistical applications the relative accuracy is of interest rather than the absolute accuracy, and these two metrics can behave quite differently. In order to characterize the dependence of the constant on the dimension of the integral, we derive the order of n needed to guarantee a vanishing relative error as the dimension of the integral increases. Formally, we wish to integrate a sequence of functions $f_p(x)$ from ${\mathbb{R}}^p\to \mathbb{R}$ with respect to a sequence of probability densities ${\varphi }_p(x)$. The Monte Carlo procedure uses the following approximation: (1) $${\int }_{D_p}{f}_p(x){\varphi }_p(x)dx\approx \sum _{i=1}^n\frac{f_p(X_i)}{n},$$(1) where X_i for $i=1.\dots ,n$ are iid samples from the probability densities ${\varphi }_p(x)$ with support $D_p\subset {\mathbb{R}}^p$. We measure the difficulty of the problem by the relationship required between p and n(p) such that when $p\to \infty$, for every $\delta >0$ (2) $$\mathbb{P}\left(\left|\frac{{\int }_{D_p}{f}_p(x){\varphi }_p(x)dx-\sum _{i=1}^{n(p)}{f}_p(X_i)/n(p)}{{\int }_{D_p}{f}_p(x){\varphi }_p(x)dx}\right|>\delta \right)\to 0,$$(2) where n(p), the sample of points being sampled, is an increasing function in p. We assume that the expectations are nonzero in order to make condition (2) well-defined. The growth rate of the function n(p) measures the complexity of the problem with respect to the dimensionality of the integral, that is, it is the sampling cost of producing an asymptotically consistent estimate in relative error. If the dimension of the integral is fixed, then consistency is guaranteed so long as $n\to \infty$. We show in the sequel that this cost can vary depending on the function being integrated, and can range from exponential to polynomial in dependency. Therefore, the performance of Monte Carlo procedures is not uniformly affected by the dimensionality of the integral.

Although considering a changing integrand may appear unusual at first, this type of scheme has been used as a theoretical device in the proof for the accuracy of quasi-Monte Carlo methods (Lemieux 2009; Owen 2019), where a carefully chosen deterministic grid is used instead of a randomly sampled grid. In this setting, uncertainty quantification is obtained through the Koksma-Hlawka inequality (Hickernell 2014), which is shown by constructing a adversarial function which changes with the number of grid points, n. This regime is also common in high-dimensional statistics, in which a statistical model is allowed to increase in complexity as the sample size increases, see Vershynin (2018) and Wainwright (2019) for an overview of high-dimensional probability and statistics, respectively. In our case, we can imagine that the function of interest lies somewhere along the sequence $f_p(x)$ and we wish to give some guidance on the required size of n to satisfy a certain accuracy requirement, depending on the placement of this function. Even though the statement in (2) is asymptotic, we use nonasymptotic concentration bounds throughout this work. Through this approach we obtain nonasymptotic $1-\alpha$ confidence intervals for the absolute error, and we derive the relationship between the sampling cost of Monte Carlo procedures in high dimensions and the structure of $f_p(x)$ and ${\varphi }_p(x)$.

We first provide some background information on Monte Carlo integration and concentration inequalities. We then examine some results for approximating volumes in high dimensions, and general integrals under structural assumptions. We then provide a comparison between asymptotic and non-asymptotic confidence intervals in Section 5. Calculations and derivations are deferred to the Appendix. We take $g(n)=\tilde{\omega }(a_n)$ to mean that $g(n)/a_n\to \infty$, and take $g(n)=O(a_n)$ to mean that there exists a constant N₀ such that for all $n>N_0$ $|g(n)/a_n|\le B$ for some B > 0. We use uppercase letters such as X or X_i to denote a vector-valued random variable. When invoking the Lipschitz property we take it to be with respect to the Euclidean norm. The Lipschitz property implies the uniform continuity of a function, and therefore bound its rate of increase uniformly across its domain. This fact will prove helpful in obtaining faster rates of concentration.

2 Background

Monte Carlo integration is justified by the application of the central limit theorem (CLT) in which we write $${\int }_{D_p}{f}_p(x){\varphi }_p(x)dx=\mathbb{E}[f_p(X)],$$where ${\varphi }_p(·)$ is a probability measure. The addition of ${\varphi }_p(·)$ makes integration on infinite domains feasible; some texts take ${\varphi }_p$ as the uniform probability measure when considering domains with finite Lebesgue measure. If we are able to sample independent and identically distributed points $X_1,\dots ,X_n$ from the density ${\varphi }_p$, we can estimate the integral through: (3) $$\mathbb{E}[f_p(X)]\approx \frac{\sum _{i=1}^n{f}_p(X_i)}{n}.$$(3)

Assuming that $E[|f_p(X)|]<\infty$, we are then guaranteed a consistent estimate as $n\to \infty$ by the strong law of large numbers.

The approximation error is typically obtained by using an empirical estimate of the variance and appealing to the CLT to produce an asymptotically valid confidence interval. However, noting that the estimate introduced in (3) is an average, quantifying the approximation error is equivalent to quantifying how rapidly an average concentrates around its expectation. Seeing the problem from this perspective allows us to use powerful concentration inequalities to bound the non-asymptotic approximation error; see Boucheron, Lugosi, and Massart (2013) for a textbook discussing concentration inequalities and their uses. In order to use these bounds however, we require additional assumptions on $f_p(x)$ and ${\varphi }_p(x)$. We introduce some particular structures which lead to rapid concentration, including averages of independent sub-Gaussian random variables, concave/convex Lipschitz functions on the hypercube and Lipschitz functions of strongly log-concave random variables.

2.1 Averages of Sub-Gaussian Variables

Sub-Gaussian random variables are a class of random variables whose tail probabilities can be dominated by those of a Gaussian random variable with mean μ and variance ${\sigma }^2$. As Gaussian distributions concentrated around their mean at an exponentially fast rate (Durrett 2019, Theorem 1.2.6), we can consequently expect the same behavior from sub-Gaussian random variables. It can be shown that for a random variable to have lighter tail than a Gaussian distribution with mean μ and variance ${\sigma }^2>0$ is mathematically equivalent to upper bounding the moment generating function of the centered random variable by that of a Gaussian with variance ${\sigma }^2$ (Wainwright 2019, Theorem 2.6). More precisely, a random variable Z is said to be sub-Gaussian with mean μ and a proxy variance ${\sigma }^2>0$ if the moment generating function of $Z-\mu$ satisfies: $$M_{Z-\mu }(t)=E\{\hspace{0.17em}\exp [t(Z-\mu )\left]\right\}\le \hspace{0.17em}\exp \hspace{0.17em}\left(\frac{{\sigma }^2{t}^2}{2}\right)\text{ for all }t\in \mathbb{R}.$$

The proxy variance ${\sigma }^2$ may not be unique in this definition, but taking ${\sigma }^2$ to be as small as possible while satisfying the above condition results in a sharper bound. Then $\overline{Z}$, the average of iid realizations $Z_1,\dots ,Z_n$ satisfies (4) $$\mathbb{P}\left(\left|\overline{Z}-E(Z)\right|>\delta \right)\le 2\hspace{0.17em}\exp \hspace{0.17em}\left(-\frac{{\delta }^{2}n}{2{\sigma }^2}\right),$$(4)

(see Wainwright 2019, sec. 2). For the purpose of verifying condition (2), the bound in (4) can be re-expressed as (5) $$\mathbb{P}\left(\left|\frac{\overline{Z}-E(Z)}{E(Z)}\right|>\delta \right)\le 2\hspace{0.17em}\exp \hspace{0.17em}\left[-\frac{{\delta }^{2}E{(Z)}^{2}n}{2{\sigma }^2}\right],$$(5) through simple algebraic manipulations. This form is more useful for studying the relative accuracy of Monte Carlo integration. We assume the expectation is nonzero or the ratio in (5) is ill-defined. As we do not know the expectation, we need to bound it from below in order to use (5) in practice; for convex functions this can sometimes be accomplished by Jensen’s inequality. One important class of sub-Gaussian random variables are bounded random variables taking value in the interval $[a,b]$, in which case the proxy variance ${\sigma }^2$ is ${(b-a)}^2/4$; this is a consequence of Hoeffding’s Lemma (Boucheron, Lugosi, and Massart 2013, Lemma 2.2).

2.2 Convex/Concave Lipschitz Functions of Random Variables on the Unit Cube

Beyond simple averages, there are concentration results for convex/concave Lipschitz functions of random variables on the hypercube. For a Lipschitz convex/concave function g(x) with Lipchitz constant L, and for a random variable $X\subset {[0,1]}^p$ whose components are independently distributed, (6) $$\mathbb{P}\left(\left|g(X)-E[g(X)]\right|>\delta \right)\le 2\hspace{0.17em}\exp \hspace{0.17em}\left(-\frac{{\delta }^2}{2L^2}\right),$$(6)

which implies (7) $$\mathbb{P}\left(\left|\frac{g(X)-E[g(X)]}{E[g(X)]}\right|>\delta \right)\le 2\hspace{0.17em}\exp \hspace{0.17em}\left\{-\frac{{\delta }^{2}E{[g(X)]}^2}{2L^2}\right\},$$(7) by Theorem 3.24 in Wainwright (2019). We will apply this bound to averages of convex/concave Lipschitz functions in Section 4, by using the fact that averages of such functions remains convex/concave Lipschitz. The surprising aspect of this bound is that it is dimension free, as it holds for all convex/concave Lipschitz continuous functions of arbitrary dimension. However, as we will see in the examples in Section 4, it is possible that the Lipschitz constant L increases with the dimension of the function.

2.3 Lipschitz Functions of Strongly Log-concave Random Variables

For some $\gamma >0$, a density $\varphi (x)=\hspace{0.17em}\exp [-\psi (x)]$ is strongly γ-log-concave if it satisfies the following: $$\lambda \psi (x)+(1-\lambda )\psi (y)-\psi [\lambda x+(1-\lambda )y]\ge \frac{\gamma }{2}\lambda (1-\lambda ){\Vert x-y\Vert }_2^2,$$for all $\lambda \in [0,1]$ and all $x,y\in {\mathbb{R}}^p$. Then for any L-Lipschitz function g and a strongly γ-log-concave random variable X: $$\mathbb{P}\left(\left|g(X)-E[g(X)]\right|>\delta \right)\le 2\hspace{0.17em}\exp \hspace{0.17em}\left(-\frac{\gamma {\delta }^2}{4L^2}\right),$$ and $$\mathbb{P}\left(\left|\frac{g(X)-E[g(X)]}{E[g(X)]}\right|>\delta \right)\le 2\hspace{0.17em}\exp \hspace{0.17em}\left\{-\frac{\gamma {\delta }^{2}E{[g(X)]}^2}{4L^2}\right\},$$ by Theorem 3.16 in Wainwright (2019). Strongly log-concave distributions can also be characterized as the product of a log-concave density with a multivariate normal distribution (Saumard and Wellner 2014). We will apply this bound to averages of iid realizations of g(X) by using properties of strongly γ-log-concave random variables in Section 4. This bound also does not depend on the dimension of the integral, although once again, the constants involved can potentially increase with p. Intuitively, the strongly log-concave density must be uni-modal and “hump-like” on the log scale, and as this is exponentiated the resulting density concentrates rapidly around its mean.

3 Estimation of High-Dimensional Volumes

Estimation of volumes is a classic application of Monte Carlo integration. Suppose we are interested in calculating the unknown volume (Lebesgue measure) of a set F which is contained within a set $F\prime$ with known finite volume. Then if we can sample X_i’s from the uniform distribution supported on $F\prime$, we can approximate the volume of F by: $$\text{Vol}(F)\approx \text{Vol}(F\prime )\sum _{i=1}^n\frac{\mathbb{I}(X_i\in F)}{n}.$$

As the random variable $\text{Vol}(F\prime )\mathbb{I}(X_i\in F)$ only takes the value of 0 or $\text{Vol}(F\prime )$, it is sub-Gaussian with ${\sigma }^2=\text{Vol}{(F\prime )}^2/4$. Noting this and using the sub-Gaussian bound in Section 2.1 gives:

Proposition 1.

Suppose that both $\text{Vol}(F)$ and $\text{Vol}(F\prime )$ are finite. Then: $$\begin{array}{c}\mathbb{P}\left(\left|\frac{\text{Vol}(F\prime )\sum _{i=1}^n\mathbb{I}(X_i\in F)}{n}-\text{Vol}(F)\right|>\delta \right)\\ \mathrm{\hspace{1em}}\mathrm{\hspace{1em}}\le 2\hspace{0.17em}\exp \hspace{0.17em}\left[-\frac{2{\delta }^{2}n}{\text{Vol}{(F\prime )}^2}\right].\end{array}$$

This implies the nonasymptotic $1-\alpha$ confidence interval $$\begin{array}{c}(\text{Vol}(F\prime )\left[\frac{\sum _{i=1}^n\mathbb{I}(X_i\in F)}{n}-\sqrt{\frac{\hspace{0.17em}\log \hspace{0.17em}(2/\alpha )}{2n}}\right],\\ \hspace{1em}\text{Vol}(F\prime )\left[\frac{\sum _{i=1}^n\mathbb{I}(X_i\in F)}{n}+\sqrt{\frac{\hspace{0.17em}\log \hspace{0.17em}(2/\alpha )}{2n}}\right])\end{array}$$for $\alpha <1$ and $$\begin{array}{c}\mathbb{P}\left(\left|\frac{\text{Vol}(F\prime )\sum _{i=1}^n\mathbb{I}(X_i\in F)/n-\text{Vol}(F)}{\text{Vol}(F)}\right|>\delta \right)\\ \mathrm{\hspace{1em}}\mathrm{\hspace{1em}}\le 2\hspace{0.17em}\exp \hspace{0.17em}\left[-\frac{2{\delta }^2\text{Vol}{(F)}^{2}n}{\text{Vol}{(F\prime )}^2}\right],\end{array}$$ where X_i are iid copies of a uniform random variable supported on $F\prime$.

The typical example provided in textbooks such as (Ross 2013, chap. 3.2), is the problem of computing the value of π or the area of a circle inscribed in a square. We extend the latter to the case of a hypersphere contained in a hypercube and showcase how the curse of dimensionality manifests itself in volume estimation.

Example 1.

Let $F=\{x\in R^p:{\Vert x\Vert }_2\le 1\}$ and let $F\prime =\{x\in R^p:{\Vert x\Vert }_{\infty }\le 1\}$, where F is a hypersphere of dimension p with radius 1, while $F\prime$ is the hypercube centered at 0 with sides of length 2. In this case we know that $\text{Vol}(F)={\pi }^{p/2}/\Gamma (p/2+1)$ and $\text{Vol}(F\prime )=2^p$, thus, $$\begin{array}{c}\mathbb{P}\left(\left|\text{Vol}(F\prime )\sum _{i=1}^n\mathbb{I}(X_i\in F)/n-\text{Vol}(F)\right|>\delta \right)\\ \mathrm{\hspace{1em}}\mathrm{\hspace{1em}}\le 2\hspace{0.17em}\exp \hspace{0.17em}\left(-\frac{{\delta }^{2}n}{2^{2p-1}}\right),\\ \mathbb{P}\left(\left|\frac{\text{Vol}(F\prime )\sum _{i=1}^n\mathbb{I}(X_i\in F)-\text{Vol}(F)}{\text{Vol}(F)}\right|>\delta \right)\\ \mathrm{\hspace{1em}}\mathrm{\hspace{1em}}\le 2\hspace{0.17em}\exp \hspace{0.17em}\left[-\frac{{\delta }^2{\pi }^{p}n}{2^{2p-1}\Gamma {(p/2+1)}^2}\right].\end{array}$$

The absolute error is decaying extremely quickly, as the volume of the unit sphere decays exponentially to 0 as the dimension increases. Therefore, the chance of hitting the sphere by sampling from the unit cube is also exponentially decaying to 0. Our estimate will essentially be 0, but this is quite close to the volume of an unit sphere in high dimensions. However, in order to satisfy our relative consistency requirement, we need the number of X_i sampled to be $n(p)=\hspace{0.17em}\exp \hspace{0.17em}\{\tilde{\omega }[p\hspace{0.17em}\log \hspace{0.17em}(p)]\}$, by Stirling’s approximation. Thus, the numerical cost of integrating this function is exponentially increasing in the dimension of the sphere.

In general this bound can be improved through a tighter control of the ratio $\text{Vol}(F)/\text{Vol}(F\prime )$. In fact, if this ratio is bounded below by a constant independent of dimension, then the relative error would be exponentially decreasing in n and we would obtain a dimension-free result. However, it is increasingly difficult to find a set $F\prime$ to cover F in high dimensions, thus, this may not be feasible in practice.

As noted by a reviewer, this example is part of the broader problem of estimating a small probability or proportion. Consider the problem of estimating the probability ${\zeta }_p:=\mathbb{P}(X\in E_p)\in (0,1]$ with the estimator ${\overline{\zeta }}_p={\sum }_{i=1}^n\mathbb{I}\{X_i\in E_p\}/n$ then: $$\mathbb{P}\left(\left|\frac{{\overline{\zeta }}_p-{\zeta }_p}{{\zeta }_p}\right|>\delta \right)\le 2\hspace{0.17em}\exp \hspace{0.17em}\left(-\frac{{\delta }^2{\zeta }_p^{2}n}{2{\sigma }^2}\right),$$by (5), where X_i are iid copies of X. The requirement for relative consistency is satisfied if ${\zeta }_p^{2}n\to \infty$, therefore, depending on the speed at which ${\zeta }_p\to 0$ an arbitrarily fast or slow rate of growth of the function n(p) is possible.

4 Expectations of High-Dimensional Functions

We consider two additional examples of Monte Carlo integration for expectations of smooth functions, going beyond the estimation of probability and proportions to demonstrate the general utility of the results presented in Section 2. The indicator function used in Section 3 is, in some sense, ill behaved, as it is not a continuous function and can change very rapidly in a small neighborhood. Considering functions with restrictions on their rate of growth through the Lipschitz property, along with other structural assumptions induces much faster rates of concentration.

4.1 Lipschitz Convex or Concave Distributions on the Unit Square

Proposition 2.

For convex or concave Lipschitz function $f_p(x):{\mathbb{R}}^p\to \mathbb{R}$ with Lipschitz constant L_p, and a random variable X whose components are independently distributed and supported on $[0,1]$, $$\mathbb{P}\left(\left|\sum _{i=1}^n{f}_p(X_i)/n-E[f_p(X)]\right|>\delta \right)\le 2\hspace{0.17em}\exp \hspace{0.17em}\left(-\frac{{\delta }^{2}n}{2L_p^2}\right).$$

This implies the nonasymptotic $1-\alpha$ confidence interval $$\left(\frac{\sum _{i=1}^n{f}_p(X_i)}{n}-\sqrt{\frac{2\hspace{0.17em}\log \hspace{0.17em}(2/\alpha )L_p^2}{n}},\frac{\sum _{i=1}^n{f}_p(X_i)}{n}+\sqrt{\frac{2\hspace{0.17em}\log \hspace{0.17em}(2/\alpha )L_p^2}{n}}\right)$$for $\alpha <1$ and $$\mathbb{P}\left(\left|\frac{\sum _{i=1}^n{f}_p(X_i)/n-E[f_p(X)]}{E[f_p(X)]}\right|>\delta \right)\le 2\hspace{0.17em}\exp \hspace{0.17em}\left\{-\frac{{\delta }^{2}nE{[f_p(X)]}^2}{2L_p^2}\right\},$$ where X_i are iid copies of X.

The proof follows from the observation that a sum of Lipschitz convex/concave functions is still Lipschitz convex/concave and the previously discussed concentration inequality. The rate of convergence for the relative accuracy only depends on the dimension through the Lipschitz constant L_p and $E[f_p(X)]$. These can increase with dimension, which places constraints on the growth rate of p relative to n, as demonstrated in Example 2.

Example 2.

Consider $${\int }_{{[0,1]}^p}{f}_p(x)dx={\int }_{{[0,1]}^p}\frac{1}{1+{\Vert x\Vert }_q}dx,$$for any q > 1. We apply Proposition 2 to show that n(p) has a polynomial dependence on p for our relative consistency condition to hold. We first note that $f_p(x)$ is $p^{\eta }$-Lipschitz, where $\eta =\max (1/q-1/2,0)$ as $$\begin{array}{c}|f_p(x)-f_p(y)|=\left|\frac{1}{1+{\Vert x\Vert }_q}-\frac{1}{1+{\Vert y\Vert }_q}\right|\\ \le {\Vert x-y\Vert }_q\le p^{\eta }{\Vert x-y\Vert }_2,\end{array}$$ since ${\Vert x\Vert }_q\le {\Vert x\Vert }_2$ for all $q\ge 2$, and ${\Vert x\Vert }_q\le p^{1/q-1/2}{\Vert x\Vert }_2$ for $1\le q<2$. The expectation of this function can be bounded from below by Jensen’s inequality: $$E\left(\frac{1}{1+{\Vert x\Vert }_q}\right)>\frac{1}{1+E({\Vert x\Vert }_q)}\ge \frac{1}{1+E({\Vert x\Vert }_1)}=\frac{2}{2+p}.$$

This results in: $$\mathbb{P}\left(\left|\frac{\sum _{j=1}^n{f}_p(X_i)/n-E[f_p(X)]}{E[f_p(X)]}\right|>\delta \right)\mathrm{}\le 2\hspace{0.17em}\exp \hspace{0.17em}\mathrm{}\left[-\frac{{\delta }^{2}n}{p^{2\eta }(2+p)}\right]\mathrm{},$$meaning that it is sufficient for $n(p)=\tilde{\omega }(p^{2\eta +1})$ to guarantee an asymptotically consistent estimator in relative accuracy.

Remark 1.

Proposition 2 applies to any function defined on a hyper-rectangle as we may shift the domain of a function to the unit cube by centering and rescaling the function, although this will change the Lipschitz constant of the function.

4.2 Strongly Log-Concave Distributions

For Lipschitz functions of strongly log-concave distributions the concentration can also be quite rapid. This allows us to consider an example of a function defined on a infinite domain and remove the convexity or concavity assumption on the function $f_p(x)$.

Proposition 3.

For a function $f_p(x):{\mathbb{R}}^p\to \mathbb{R}$ with Lipschitz constant L_p, and a strongly γ-log-concave random variable X, $$\mathbb{P}\left(\left|\sum _{i=1}^n{f}_p(X_i)/n-E[f_p(X)]\right|>\delta \right)\le 2\hspace{0.17em}\exp \hspace{0.17em}\left(-\frac{\gamma {\delta }^{2}n}{4L_p^2}\right).$$

This implies the nonasymptotic $1-\alpha$ confidence interval $$\left(\frac{\sum _{i=1}^n{f}_p(X_i)}{n}-\sqrt{\frac{4\hspace{0.17em}\log \hspace{0.17em}(2/\alpha )L_p^2}{\gamma n}},\frac{\sum _{i=1}^n{f}_p(X_i)}{n}+\sqrt{\frac{4\hspace{0.17em}\log \hspace{0.17em}(2/\alpha )L_p^2}{\gamma n}}\right)$$for $\alpha <1$ and $$\mathbb{P}\left(\left|\frac{\sum _{i=1}^n{f}_p(X_i)/n-E[f_p(X)]}{E[f_p(X)]}\right|>\delta \right)\le 2\hspace{0.17em}\exp \hspace{0.17em}\left\{-\frac{\gamma {\delta }^{2}nE{[f_p(X)]}^2}{4L_p^2}\right\},$$ where X_i are iid copies of X.

This is proven by noting that an average of Lipschitz function is Lipschitz and the fact that the joint distribution of independent strongly γ-log-concave distribution is still strongly γ-log-concave, as shown in the Appendix. One particularly important member of strongly log-concave distributions is the multivariate normal distribution with nonsingular covariance matrix Σ, which is strongly γ-log-concave, with γ equal to the minimal eigenvalue of ${\Sigma }^{-1}$.

Example 3.

Consider the following expectation where $X\in {\mathbb{R}}^p$ follows a multivariate normal distribution with 0 mean and a covariance matrix Σ such that ${\sigma }_{\min }({\Sigma }^{-1})=\gamma >0$. This condition on the minimum eigenvalue implies that the distribution of X is strongly γ-log-concave. We wish to estimate: $$\mathbb{E}\left[\arctan \left(1+{\Vert X\Vert }_1\right)\right].$$

This function is $p^{1/2}$-Lipschitz due to the fact that: $$\arctan (a)\text{-}\arctan (b)\le \frac{a-b}{1+ab},$$and the fact that $\Vert x{\Vert }_1-\Vert y{\Vert }_1\le \Vert x-y{\Vert }_1$. Finally we may bound the expectation of this function from below by $\arctan (1)$, thus, $$\mathbb{P}\left(\left|\frac{\sum _{i=1}^n{f}_p(X_i)/n-E[f_p(X)]}{E[f_p(X)]}\right|>\delta \right)\le 2\hspace{0.17em}\exp \hspace{0.17em}\left[-\frac{\gamma {\delta }^{2}n\arctan {(1)}^2}{4p}\right],$$ meaning that $n(p)=\tilde{\omega }(p)$ is sufficient to guarantee a relative consistent estimate.

5 Nonasymptotic and Asymptotic Confidence Intervals

We briefly discuss and compare the non-asymptotic confidence intervals presented in Sections 3 and 4 to asymptotic and exact confidence intervals in a simple example to better understand their coverage properties. The nonasymptotic confidence intervals considered in the previous sections will tend to be conservative for most distributions. This is necessary as the probability bounds used are applicable to a wide class of distributions and the equality in the bound cannot hold for all values of $t\in \mathbb{R}$, thus, they must in general overstate the probability of the event.

Example 4.

We observed a sample from a binomial distribution with parameters (k, p) where k is the number of coin tosses and p is the probability of heads. We treat k as known and consider the problem of providing a 90% and 95% confidence interval for p. We consider three possible confidence intervals, the Clopper-Pearson exact confidence interval, the Bayesian credible interval obtained with the $\text{Beta}(1/2,1/2)$ prior (which is the Jeffreys’ prior) and the nonasymptotic confidence interval from Proposition 1. The Bayesian credible interval is used instead of the usual normal approximation, as for certain values of p we are likely to obtain an observation of 0 or k, which results in an invalid approximation. This credible interval also coincides with the standard normal approximation as $k\to \infty$ by the Bernstein von-Mises theorem (Van der Vaart 2000, chap. 10), and is therefore asymptotically valid. Exact coverage probabilities for 90% and 95% nominal confidence intervals are given in Figures 1 and 2. We see that the proposed nonasymptotic interval tends to be more conservative than the Clopper-Pearson interval. Although, we do note the Bayesian credible interval suffers poor coverage for certain values of the parameter p even if the sample size is relatively high.

Fig. 1 Exact coverage probabilities for a nominal 90% confidence interval, for k = 10 on the left and k = 100 on the right. Non-asymptotic coverage probabilities are given by triangles, Clopper-Pearson exact confidence interval coverage probabilities are given by circles and Bayesian credible interval coverage probabilities are given by squares.

Fig. 2 Exact coverage probabilities for a nominal 95% confidence interval, for k = 10 on the left and k = 100 on the right. Nonasymptotic coverage probabilities are given by triangles, Clopper-Pearson exact confidence interval coverage probabilities are given by circles and Bayesian credible interval coverage probabilities are given by squares.

In general, if sufficient knowledge of the distribution is available to the practitioner, it is possible to obtain better alternatives to the non-asymptotic intervals by performing a bespoke analysis. In the binomial example, exact intervals are available as we have the analytical form of the probability mass function. However, it is not always the case that such knowledge available to us. Although, we do note that if an exact interval is unavailable, the asymptotic confidence intervals may produce invalid anti-conservative statements while the nonasymptotic intervals provides valid conservative statements.

6 Discussion

We have shown that the applicability of Monte Carlo integration varies substantially depending on the specific combination of the function being integrated and the measure it is integrated against, and we relate the approximation error to the rate of concentration of a sum of iid random variables. We note that the procedure can be catastrophically bad for approximating volumes, but for integrating specific functions it may perform well. The analysis can potentially be extended to sub-exponential random variables, or random variables defined on manifolds by exploiting the numerous available concentration inequalities. Another important extension is to consider non iid random variables. This generalization could help in providing non-asymptotic bounds on the performance of Monte Carlo Maximum Likelihood Estimates for exponential graphical models as introduced in Geyer and Thompson (1992) and general Monte Carlo maximum likelihood estimation approaches. This can potentially be done through the bounds of strongly log-concave densities presented in Section 4.2, however, we would need to ensure that the sampled distribution is strongly log-convex, which is a nontrivial assumption to verify.

There are also other variants of Monte Carlo which are used in high dimensions. While quasi-Monte Carlo methods has a absolute error rate of $O[\log \hspace{0.17em}{(n)}^p/n]$, these methods perform better than the $\hspace{0.17em}\log \hspace{0.17em}{(n)}^p$ factor in the error might suggest. For example, it has been observed that quasi-Monte Carlo performs well in some financial applications in very high dimensions (Paskov and Traub 1996). This may be due to the fact that this error bound is obtained through a worst case scenario analysis and the true dependence on dimension may be much smaller. Determining the effective dimension of the problem is still an open question, see Owen and Pan (2021).

Supplementary Materials

The supplementary materials contains the code used to produce Figures 1 and 2.

Acknowledgments

We would like to thank two anonymous reviewers and the associate editor for their comments which lead to improvements in the quality and exposition of the article. We also thank Michaël Lanlancette and Heather S. Battey for comments and feedback on a draft version of this article.

Disclosure Statement

The author reports there are no competing interests to declare.

Appendix A.

Appendix

A.1 Properties of Lipschitz and Log-Concave Distributions

Lemma 1.

• If $f_i(x_i):{\mathbb{R}}^p\to \mathbb{R}$ are L_p-Lipschitz functions with respect to the Euclidean norm for $i=1,\dots ,n$, then $g(x_1,\dots ,x_n)={\sum }_{i=1}^n{f}_i(x_i)/n$ is a $L_p/n^{1/2}$-Lipschitz function.

• If $f_i(x_i):{\mathbb{R}}^p\to \mathbb{R}$ are L_p-Lipschitz convex/concave functions with respect to the Euclidean norm for $i=1,\dots ,n$, then $g(x_1,\dots ,x_n)={\sum }_{i=1}^n{f}_i(x_i)/n$ is a $L_p/n^{1/2}$-Lipschitz convex/concave function.

Proof.

Let $x_i,y_i\in {\mathbb{R}}^p$, then $$\begin{array}{c}\left|\frac{1}{n}\sum _{i=1}^n{f}_i(x_i)-\frac{1}{n}\sum _{i=1}^n{f}_i(y_i)\right|\le \frac{1}{n}\sum _{i=1}^n|f_i(x_i)-f_i(y_i)|\\ \le L_p\frac{\sum _{i=1}^n{\Vert x_i-y_i\Vert }_2}{n}\\ =L_p\frac{\sum _{i=1}^n{\left[{\Vert x_i-y_i\Vert }_2^2\right]}^{1/2}}{n}\\ \le \frac{L_p}{n^{1/2}}{\Vert (x_1,\dots ,x_n)-(y_1,\dots ,y_n)\Vert }_2,\end{array}$$where the last inequality follows from Jensen’s inequality for the concave function $\rho (x)=x^{1/2}$. The second statement follows from the first by noting that sums of convex/concave functions remain convex/concave. □

This basic property of strongly log-concave distribution is shown in Saumard and Wellner (2014), but due to notational differences we provide a proof for the convenience of the reader.

Lemma 2.

Let X_i be independent copies of a γ strongly log-concave distribution for $i=1,\dots ,n$, then the joint distribution $(X_1,\dots ,X_n)$ is also a γ strongly log-concave distribution.

Proof.

Let X_i be distributed according to ${\varphi }_i(x_i)=\hspace{0.17em}\exp [-{\psi }_i(x_i)]$ Let $y_i,z_i\in {\mathbb{R}}^p$, for $i=1,\dots ,n$ then by the log-concavity of each x_i, for each $i=1,\dots ,n$ and all $\lambda \in [0,1]$ $$\begin{array}{c}\lambda {\psi }_i(y_i)+(1-\lambda ){\psi }_i(z_i)-{\psi }_i[\lambda y_i+(1-\lambda )z_i]\\ \ge \frac{\gamma }{2}\lambda (1-\lambda ){\Vert y_i-z_i\Vert }_2^2,\end{array}$$which implies $$\begin{array}{c}\lambda \sum _{i=1}^n{\psi }_i(y_i)+(1-\lambda )\sum _{i=1}^n{\psi }_i(z_i)-\sum _{i=1}^n{\psi }_i[\lambda y_i+(1-\lambda )z_i]\\ \ge \frac{\gamma }{2}\lambda (1-\lambda )\sum _{i=1}^n{\Vert y_i-z_i\Vert }_2^2\\ \ge \frac{\gamma }{2}\lambda (1-\lambda ){\Vert (y_1,\dots ,y_n)-(z_1,\dots ,z_n)\Vert }_2^2,\end{array}$$ where the final inequality follows from the triangle inequality. □

A.2 Proof of Propositions 2 and 3

A.2.1 Proof of Proposition 2

The first and third part of the proposition follow from choosing the function $g(X_1,\dots ,X_n)={\sum }_{i=1}^{n}f(X_i)/n$ in bound (5), and noting by Lemma 1 that g is convex/concave with Lipschitz constant $L_p/n^{1/2}$. The confidence interval is obtained by setting: $$\alpha =2\hspace{0.17em}\exp \hspace{0.17em}\left(-\frac{{\delta }^{2}n}{2L_p^2}\right),$$and solving for δ.

A.2.2 Proof of Proposition 3

The proof for the first and third statements is almost the same as that of Proposition 2. We let $g(X_1,\dots ,X_n)={\sum }_{i=1}^{n}f(X_i)/n$ and note that this function is $L_p/n^{1/2}$-Lipschitz and combine this with the fact that the joint distribution of $(X_1,\dots ,X_n)$ is also strongly γ-log-concave (as shown in Lemma 2) gives the desired result. The confidence interval is obtained by setting: $$\alpha =2\hspace{0.17em}\exp \hspace{0.17em}\left(-\frac{\gamma {\delta }^{2}n}{4L_p^2}\right),$$and solving for δ.