The bootstrap method for Monte Carlo integration inference

1 Introduction

In the real world applications, the standard or classical mathematical methods have been widely used to compute the exact integral of a function g defined on an interval $[a,b]$ . This can be achieved if the situation is simple; however, if the situation is complicated, the computation is quite hard or impossible in some cases. Therefore, the approximation methods can be good choices to compute the integrals for the complicated functions with minor errors. One of these approximation methods is the Monte Carlo integration method, which is described in many references, see e.g. (Kalos and Whitlock, 2009; Rizzo, 2019; Yang, 2014) for more details.

This approximation method for integrals is built with respect to the sampling distribution. It is crucial to choose a suitable distribution, which should be close to the function g, to generate random observations. The suitable distribution can lead to have a better approximation for the integral with small variance. Many references in the literature discussed the importance sampling or variance reduction, see e.g. (Hammersley and Morton, 1956; Oh and Berger, 1992; Rizzo, 2019; Rubinstein and Marcus, 1985; Tokdar and Kass, 2010; Van Dijk and Kloek, 1983) for detailed presentations.

Surely, the importance sampling is crucial for the Monte Carlo integration to derive a well approximation with small variance, but the importance sampling requires a long computational time and some information about the shape of g. These requirements motivate to use the bootstrap method with the Monte Carlo integration. Based on the bootstrap distribution, it is possible to derive well approximations for integrals with small standard errors. Moreover, the percentile and Normal confidence intervals of any integral can be derived with high accuracy.

This paper is organized as follows: Section 2 presents the Monte Carlo integration for one dimension along with some descriptions. In Section 3, the bootstrap method is described with explanations of computing the standard error and deriving the percentile and Normal confidence intervals for the integral of interest $\theta$ . Section 4 presents the performance of bootstrap with Monte Carlo integration. The last section presents some concluding remarks.

2 Monte Carlo integration

The Monte Carlo method is a well-known concept used to derive the approximate integrals for complicated functions defined on certain domains. This method is built based on the probability theorem (Chung and Zhong, 2001; Jaynes, 2003; Loeve, 2017), where the integral of a function g can be computed by taking the expectation of function $g,E\left\{g\right\}$ . For univariate data, suppose that X is a random variable following a probability density function f. The expectation of $g\left(X\right)$ can be written mathematically as follows

Let $\theta ={\int }_a^{b}g\left(x\right)\mathit{dx}$ and the random variables $X_1,X_2,\dots ,X_n$ are independent and identically distributed from the Uniform distribution with parameters a and b. Further suppose that $x_1,x_2,\dots ,x_n$ are the corresponding observations to the random variables $X_1,X_2,\dots ,X_n$ . The estimated result of $\theta$ can be found by

1. Generate n observations from the Uniform distribution with parameters a and b.

2. Find $g\left(x_i\right)$ for all $i=1,2,\dots ,n$ .

3. Compute $\widehat{\theta }=\frac{b-a}{n}{\sum }_{i=1}^{n}g\left(x_i\right)$ .

The Uniform distribution can be replaced by any known distribution to generate the observations in the first step, but we choose the Uniform distribution here for simplicity in application.

3 The bootstrap method

The bootstrap method is a resampling concept proposed to measure the accuracy of a statistical estimate and to make inferences about unknown population parameters, e.g. mean, median, variance and confidence interval. In the literature, it has been widely used due to its simplicity to apply and efficiency to provide well estimates. With high accuracy, the bootstrap distribution for any statistic of interest can mimic the sampling distribution, where the sampling distribution is not always easy to obtain in real applications.

To implement the bootstrap method, parametric and nonparametric models are used to create multiple bootstrap samples, then the statistic of interest is computed based on each bootstrap sample. The empirical distribution of the results can be used as a proxy distribution for the sampling distribution and this allows making inferences about the statistic of interest. For more detailed presentation, it is beneficial to see the book of “An Introduction to The Bootstrap” by Efron and Tibshirani (1993) and the book of “Bootstrap Methods and Their Application” by Davison and Hinkley (1997).

To use the bootstrap method with Monte Carlo integration, we first need to introduce some notations. Suppose that the random quantities $X_1,X_2,\dots ,X_n$ are independent and identically distributed following the probability distribution f and supported on $[a,b]$ . Let $x_1,x_2,\dots ,x_n$ be the observations corresponding to these random quantities. Furthermore, let the integral of interest be $\theta ={\int }_a^{b}g\left(x\right)\mathit{dx}$ . Now, it is easy to present the algorithm of the bootstrap method with Monte Carlo integration through the following steps.

1. Generate n observations from the probability distribution f.

2. Find $g\left(x_i\right)$ for all $i=1,2,\dots ,n$ .

3. Compute ${\widehat{\theta }}^{\ast }=\frac{b-a}{n}{\sum }_{i=1}^{n}g\left(x_i\right)$ .

4. Perform steps (1), (2) and (3) B times; this leads to ${\widehat{\theta }}^{*1},{\widehat{\theta }}^{*2},\dots ,{\widehat{\theta }}^{*B}$ .

For better performance, the value of B is suggested to be large; e.g.  $B=1000$ (Al Luhayb, 2021; Al Luhayb et al., 2023; Efron, 1979, 1981). To derive a bootstrap estimate of $\theta$ , we compute the average of ${\widehat{\theta }}^{*1},{\widehat{\theta }}^{*2},\dots ,{\widehat{\theta }}^{*B}$ , and to provide a bootstrap standard error estimate for $\overline{{\widehat{\theta }}_{\mathit{boot}}^{\ast }}=\frac{1}{B}{\sum }_{j=1}^B{\widehat{\theta }}^{*j}$ , we compute the standard deviation of ${\widehat{\theta }}^{*1},{\widehat{\theta }}^{*2},\dots ,{\widehat{\theta }}^{*B}$ by Efron and Tibshirani (1986)

For the $(1-\alpha )\%$ percentile confidence interval, we order the values ${\widehat{\theta }}^{*1},{\widehat{\theta }}^{*2},\dots ,{\widehat{\theta }}^{*B}$ from least to largest, then take the $\left(\frac{\alpha }{2}\right)$ th and $(1-\frac{\alpha }{2})$ th ordered values, where the former is the lower bound and the latest is the upper bound. This can be written as follows

For the $(1-\alpha )\%$ Normal confidence interval, it is needed to compute $\overline{{\widehat{\theta }}_{\mathit{boot}}^{\ast }}$ and $\mathit{SE}\left(\overline{{\widehat{\theta }}_{\mathit{boot}}^{\ast }}\right)$ , then we use the following equation (Hazra, 2017)

4 Simulation studies

In this section, we present different bounded integrals needed to be computed analytically based on the bootstrap method with Monte Carlo integration. We choose different domains and functions as shown in Fig. 1. Also, we determine the Uniform and Normal distributions, restricted to the domain of the integral of interest, to generate observations for the bootstrap samples. To compute the bootstrap estimate of integral $\theta$ with the standard error and the percentile and Normal confidence intervals, we set B equal to 1000 and the sample size of each bootstrap sample n equal to 1000. This means that the iteration of bootstrap is 1000 as $B=1000$ . Table 1 presents different bounded integrals with their analytical results, and these integrals will be estimated based on the bootstrap method with Monte Carlo integration to make comparisons with the true results. By this strategy, we can make investigations on the performance of our method.

Close

Fig. 1

The shape of each function $g\left(x\right)$ presented in Table 1.

Export to PPT

Tables 2–5 present the bootstrap estimates $\widehat{\theta }$ along with the standard errors of $\widehat{\theta }$ and the $90\%$ percentile and Normal confidence intervals of $\theta$ for all examples presented in Table 1. It is obvious that the bootstrap estimates are nearly identical to the exact results. This is the power of using the bootstrap method with Monte Carlo integration, which is more needed for complicated integrals that is impossible to be computed theoretically. Through the bootstrap procedure, we use Uniform and Normal distributions with different parameters, but restricted to the integral’s bounds, and this leads to have nearly identical standard errors. From this observation, we can assure that the bootstrap distribution can be a good proxy distribution for the integral $\theta$ ; there is no need to estimate the shape of $g\left(x\right)$ for importance sampling, or variance reduction. This helps to conserve time to running codes and it is possible to do analysis with less information. Also, based on the bootstrap method, the $90\%$ percentile and Normal confidence intervals can be easily derived. From Tables 2–5, we note that the $90\%$ percentile confidence intervals are all wider than the $90\%$ Normal confidence intervals.

5 Concluding remarks

In this paper, we illustrated the bootstrap method with Monte Carlo integration for one dimension, and the method was used through multiple examples. From the examples, it can be concluded that the bootstrap method with Monte Carlo integration is a good approach to compute well approximate integrals for different functions in different domains with small variances, regardless to the distribution being used to generate observations for the bootstrap samples. Based on the bootstrap estimates, the $(1-\alpha )\%$ percentile and Normal confidence intervals for $\theta$ can be derived with high accuracy. Our method can be beneficial for integrals that may be difficult or impossible to compute. Also, the method is easy to implement and straightforward, which is only relying on sampling from a known distribution, then taking the evaluation. This is repeated multiple times and this should be large, e.g.  $B=1000$ . To put the method into a practical use, we included the R codes in the appendix. To run the codes, it requires about ten seconds, which is nothing in real applications and this is one of the advantages.

As a future research, the method will be generalized for multiple dimensions with more complicated functions and domains. To achieve this generalization, we may use the copula concept, which is able to take the dependence structure between the variables into account, to generate observations for the bootstrap samples. For more detailed presentations about the copula concept, it is advised to see Coolen-Maturi et al. (2016); Muhammad (2016); Muhammad et al. (2016); Muhammad et al. (2018); Sklar (1959).